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A simple model of the non-equilibrium dissociation 
of a gas in Couette and boundary-layer flows 


By JAMES E. BROADWELL 
Santa Monica Division, Douglas Aircraft Company, California* 


(Received 11 September 1957 and in revised form 10 February 1958) 


SUMMARY 

At sufficiently high temperature the oxygen and nitrogen 
molecules in air dissociate into atoms. Energy is required for 
this decomposition and, furthermore, if the temperature is not 
uniform, concentration gradients are formed and energy is 
transferred by the consequent interdiffusion of atoms and 
molecules. In this paper, a simple model of a gas is formulated 
to illustrate the effect of these processes on heat transfer in three 
situations: (1) ina layer of gas at rest between two walls at different 
temperature, (2) in Couette flow, and (3) in a laminar boundary 
layer on a flat plate. 

In the model, the law for the rate of reaction (dissociation and 
recombination) is of an especially simple form with the speed of 
reaction characterized by a single parameter. When the parameter 
becomes large, the gas approaches chemical equilibrium; at the 
other extreme, no reaction occurs within the gas. 

Two chemical conditions of the walls are considered, one being 
that the walls are catalytic (the surface reaction rate is assumed 
sufficiently high to hold the gas at the wall in chemical equilibrium), 
and the other that the walls have no effect on the reaction, 1.e. they 
are non-catalytic. 

The most important of the simplifications made are: (a) The 
reaction rate law is put in a form in which the equilibrium 
concentration of atoms varies linearly with temperature. ‘Thus, 
there is only one temperature at which the gas is undissociated 
instead of the actual range of temperature. (6) In problems (1) 
and (2), », k, and pD are taken to be constant. (c) In problem (3) 
the Lewis number is assumed to be unity. (d) Only binary 
mixtures of atoms and molecules are considered (so that ionization 
and more complex dissociation processes are not covered). 
(e) Coupling of irreversible flows, such as that causing thermal 
diffusion, is neglected. (f/f) Only problems in which the pressure 
is constant are considered. 


* The work in this paper was performed while the author was on leave from the 
Department of Aeronautical Engineering, University of Michigan. 


F.M. H 
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1. INTRODUCTION 

When the temperature of air is raised sufficiently, for instance, by 
being brought to relative rest in the boundary layer on a body moving at 
a very high speed, some of the molecules of the constituent gases dissociate 
into atoms. Energy is ‘absorbed’ in the dissociation and is transported 
by the interdiffusing atoms and molecules. Both processes affect the heat 
transfer through the mixture. This paper describes three instances of 
heat transfer in a simple model of a dissociating (and recombining) gas. 
These are: (1) ‘The heat transfer through a layer of the gas at rest between 
two walls at different temperature. (2) Couette flow. (3) Boundary-layer 
flow. 

Leipmann & Bleviss (1956) and Liepmann & Roshko (1957) have 
considered the Couette flow of a dissociated gas which is in chemical 
equilibrium. ‘Their work suggested the present study in which the 
assumption of equilibrium is dropped but in which other simplifications 
are introduced. 

A paper by Hirschfelder (1957) dealing with heat transfer in a chemically 
reacting mixture appeared after the present work was essentially complete. 
There the exact equations governing the heat transfer through a layer of 
a reacting gas are discussed and an iterative method of solution, applicable 
when the reaction rates are high, is developed. The similarity between the 
character of these solutions and those of the present model, when a reaction 
rate parameter is large, are pointed out below. 


1.0| 


a, 








Tt 


Figure 1. Equilibrium dependence of atom concentration on temperature. 


To see how the heat transfer process is modified by the dissociation 
of the fluid, consider the forward and reverse reactions O, = 20 taking 
place in a closed vessel. These reactions establish an equilibrium 
composition that is a function of pressure and temperature. Thus, if the 
gas is heated slowly and uniformly at constant pressure, the mass fraction, 
a,, of atoms varies as shown in figure 1 (see, for example, Lighthill (1957)). 

As the atoms appear, the curve of enthalpy vs temperature becomes 
steeper; this is a consequence of the energy required to decompose the 


molecule. 
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On the other hand, if the gas is not heated uniformly, spacewise variations 
in composition as well as in temperature appear. ‘The temperature gradients 
establish heat fluxes and the concentration gradients give rise to mass flows 
and diffusion. ‘The diffusing components transport their enthalpy with them, 
and this transport constitutes an energy flow that may be large compared 
with that caused by the thermal gradients. 

In general then, dissociation (and recombination) directly alters heat 
transfer processes: firstly by creating ‘sources’ and ‘sinks’ of energy, and 
secondly by providing, through diffusion currents, another energy transport 
mechanism. 

There are, of course, other indirect (but important) consequences, such 
as the change in viscosity, specific heats, gas constant, and so forth. Let 
us consider again a gas in an enclosure. Suppose the temperature is 
suddenly and uniformly increased, then the rate of production of atoms 
exceeds the rate of disappearance and the mixture approaches the new 
equilibrium at a finite rate which depends upon the pressure, the temperature, 
and the concentration of atoms. Now notice that a gas can be kept 
permanently out of equilibrium by diffusion, even when it is stationary 
and the heat transfer processes are steady. It is, therefore, of interest 
and it is one purpose of this paper to study how departure from equilibrium 
affects the heat transfer mechanism in a dissociating gas. 

In this paper only the dissociation of a single diatomic gas into a binary 
mixture of atoms and molecules is treated. ‘The ionization and dissociation 
of more complex gases and mixtures of gases is not considered. 

In § 2 the equations governing the flow of a reacting gas are given. The 
simplifications leading to the model gas are discussed in § 3 and its behaviour 
in the three situations described above is presented in the final sections. 


2. GENERAL EQUATIONS 


2.1. Conservation equations 

The equations governing the flow of a mixture of interdiffusing and 
chemically reacting gases are derived in Hirschfelder, Curtiss & Bird (1954). 
Before these equations can be stated, however, it is necessary to define the 
several velocities that arise in such a flow. 

Consider a mixture of gases moving relative to some fixed coordinate 
system. A small plane, perpendicular to the flow of constituent 7, and 
moving so that, on the average, no molecules of 7 cross it, moves with the 
velocity of 7, y;. ‘The mass flow of 7 per unit area relative to the fixed 
coordinate system is m;m, y;, where n is the number of moles per unit volume 
and m is the molecular weight. The total mass flow per unit area is 
n,m; y;. Thus, with the mass density p given by p= }n;m,, the mass 
average velocity V is defined by V=(l/p)Sn,m,y,. Finally, the 
diffusion velocity of the ith component, V,, is the velocity of that component 
relative to the mass average velocity; thus V,; = y;—V. 

H2 
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In the following, the conservation equations, taken from Hirschfelder 
et al. (1954), are stated first for a mixture composed of an arbitrary number 
of components and then for a mixture containing but two. In the problems 
treated later, these two components will be the atoms and molecules of a 
dissociated diatomic gas. 

Only very simplified forms of these equations are needed in the first 
two problems treated in this paper. Nevertheless, we give the complete 
equations here in order to describe more clearly the physical processes 
involved and because they are needed in the boundary layer problem of § 6. 


Conservation of mass 
The equation for the conservation of mass of each species can be written 
on,/ot + V.n(V+V;) = K;, (1) 
where K; is the mole rate of production of the 7th component per unit 
volume by chemical reaction. 

Note first that the multiplication of this equation by m,, the summation 
of the resulting equation over all 7, and use of the overall mass conservation 
condition > m;K; = 0, lead to the usual continuity equation 

Cp/ot + V. (pV) = 0. (2) 

The variables V; and A, in equation (1) are determined by the 
thermodynamic properties and their gradients. Consider first the diffusion 
velocity. In general, diffusion may be caused by gradients in the pressure 
and temperature as well as in concentration. In many important situations 
the last-named cause is dominant, and it is the only one considered in the 
following discussion, i.e. cross-coupling of the fluxes is neglected. If, in 
addition to this restriction, the mixture consists of only two components, 
then (see Hirschfelder et al. (1954, p. 516)) 

V, (D/x;)Vz;, $= 1, Z, (3) 
where D is the binary diffusion coefficient and z; the mass fraction of the 
ith component defined by 


n; mM; p; 
x $< =~. 
N,m,+n,m, ~ p 
We turn next to the variable A;, the net rate of production of species 7 
by chemical reaction. Such a reaction is described by an equation of the 
form 
Bid, + Botot... = yd + rvedot... 
J; denotes the ith chemical component, and § and v the stoichiometric 
coethcients of reactants and products, respectively. The net rate of 
reaction, r, can be put in the form (see Hirschfelder et al. 1954, p. 748) 


r=K,(T)nin>...—K,(T)nine..., 


in which A, and K, are the rate constants for the forward and reverse 


j 


reaction. ‘Then 
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Also, introducing the mole fraction, X; = n,;/n and with n = p/(RT), we 
may write 


Bi 2; 
_ KD fr) XpXB...-K,(7)( £5) XnX... 


As a specific example of a reaction rate law, consider that for the 
decomposition of oxygen O, = 20. Hirschfelder (1957) states that in 
this reaction the net rate of production of atoms, Kj, is given by 

K, = —n®K,[X?—{(1—X,)/pjexp(1-58 — 60000/7)], 
where 
K, = 1:2 x 1016(7/300)*/? cm®/mole? sec, 
and p is pressure in atmospheres, 7 the absolute temperature in degrees 
Kelvin, and Y, denotes the mole fraction of atoms. 

It should be noted that the equation formed by setting A; equal to 
zero yields the equilibrium composition of the mixture as a function of the 
pressure and temperature. 

The expression for the reaction rate A; and the above equation for the 
diffusion velocity V; may now be put in equation (1) to form the mass 
conservation expression for the ith component. When we have two 
components we need, in addition to the overall continuity equation (2), 
a conservation equation for only one of the two, for «, say, 

p(ca,/et + V.Va,)— V. (pD Va,) = m, Kj. (4) 
This equation comes, with some algebraic manipulation, from equation (1) 
with the use of equations (2) and (3). ‘The appropriate expressions for A, 
will be inserted later. 
Conservation of momentum 

In terms of the mass average velocity, V, with components w,, the 
equations of motion are the usual ones* : 


Du 2 2 a Cu, 0 Ou Ou, ‘ 
p -_=- = —=>=—(u~—]+ —|Bv 7 = + =} |, (5) 
Dt Om 30%, OX, OX, OXp OX, 


x 
in which , is the viscosity of the mixture and 
D/Dt = 0/et +u, 0/dx,. 








Conservation of energy 
The energy equation is 
Dh Dp | 
"Dt Dt 


in which ® is the dissipation function, given by 


eu,\(Cu, , OU, 2 ja. 
© ~ ol Se )(ae + ae) — eae) 
0X2] \C X 5 OX, < OX, 


* The equation given by Penner (1955) and derived by von Karman contains 
additional stress terms arising from diffusion. In general, these would be small 
compared to the usual stress terms and need not, especially in view of the later 


V.(R VT — > p;h, V,)+9, (6) 








approximations, be included here. 
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and h is the enthalpy per unit mass of mixture. The energy carried by 
diffusion is }p;h;V, in which 4h; is the enthalpy per unit mass of 
component 7. The symbol k denotes the thermal conductivity of the 
mixture with the components treated as non-reacting gases. 


2.2. Thermodynamic properties of the mixture 


The equation of state of a mixture of perfect gases is 
p a (> R; *) 7; 


which for a binary mixture becomes 


R ms ni 
=p—|]1+(—-1 T. 
P al 6 Ja | 


R, is the gas constant for the ith component. 
The enthalpy per unit mass of mixture is given by h = }2;h,, in 

which 

-T 

h,= | (C,)dT+h), 

40 
where A? is the enthalpy of formation extrapolated to absolute zero and C,, 
the specific heat at constant pressure. 


2.3. Boundary conditions 


The boundary conditions that are to be applied to the temperature 
and mass average velocity are the usual ones and require no special 
discussion. ‘They will be stated later for each problem. However, the 
conditions to be satisfied by the concentrations at solid boundaries may 
not be equally familiar. They are discussed briefly here. 

If the wall has a catalytic effect on a reaction, then there may be a net 
rate of production of the species involved, per unit area of the wall. In 
steady flow, this production of each species must be balanced by the flux 
normal to, and at, the wall. In the following, two simple extremes are 
considered: (1) the wall is not catalytic, and (2) the wall reaction rate is 
such as to hold the concentrations there in chemical equilibrium at the 
wall temperature. 

When the wall is non-catalytic, the flux of each component normal to 
the wall must vanish at the wall. Hence, if, as in the problems considered 
in this paper, the gradients are normal to the wall (or in the boundary-layer 
case are assumed to be), then 


n; m, V, wall = — pDVx, wall = 0. 


3. IDEALIZED GAS AND SIMPLIFIED EQUATIONS 


In this section we state the properties of a simplified model of a gas that 
approximates, at least qualitatively, the behaviour of the mixture of atoms 
and molecules formed when a gas, such as oxygen, is subjected to a 
temperature gradient at temperatures high enough to cause dissociation. 
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First, the two components are assumed to be perfect gases with constant 
specific heats. Hence 
h, = (C,), T +h, hy = (Cy)eT +h, 
where h° is the enthalpy of formation and the subscripts 1 and 2 denote 
atoms and molecules, respectively. ‘Then 
h= ah, +h, 
or, since %,+% = l, 
h = (C,)o ts a a [(C,)s a (C,,)o] T+ Ah x,, 
where Ah® = hi—hj. That is, Ah® is the enthalpy of dissociation at 
absolute zero and A} is taken to be zero. Now {(C,),—(C,)s}/(C,)e2 1s 
small; it is of the order of 10-! for oxygen. Further, [(C,),—(C,)o]T 
is small compared with Ah° for temperatures not greatly exceeding those 
for complete dissociation. Hence, we neglect the second term in comparison 
with the first and third and drop the subscript on C,, to obtain 
h=C,T+Ah°x,. 











Figure 2. Reaction rate law. 


Consider next the reaction rate law, an example of which was given 
in §2. When a single gas dissociates or recombines and the pressure is 
constant, then the net rate of formation of atoms, K,, is a function of the 
temperature and the weight fraction of atoms, i.e. K, = K,(T,«,). The 
shape of this function depends, of course, on the particular reaction but 
in general is as sketched in figure 2. 
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‘The approximation to the reaction rate law that is to be used in this 
paper consists of the first-order terms in a Taylor series expansion of K, 
about a point O on the equilibrium line 


OR ve ws fom 7 
K, = (SF) (=~ 1) +(5 ) (%, — %19)- (7) 
C 0 0%, Jo 


(In the boundary layer problems, §6, this approximation is slightly 
modified.) This planar approximation for the surface A, = A,(7,~,) 
cannot, of course, be accurate over the entire range of the variables. It 
requires, for instance, a linear variation of the equilibrium concentration 
of atoms with temperature, a result that is obviously incorrect near «, = 0 
and x, = 1-0, and which, of course, makes no sense when the temperature 
is either below that at which dissociation begins or exceeds that for complete 
dissociation. Furthermore, in the case of oxygen, while it is true that the 
equilibrium concentration of atoms is approximately linear with temperature 
for x, not near zero or unity, the slope, 0K,/¢cx,, of the surface varies greatly 
along the equilibrium line. A rate law of the form in equation (7) does, 
nevertheless, describe qualitatively the dissociation (and recombination) 
at a finite rate of a gas in a range of temperature between the onset and 





completion of dissociation. 

The stationary gas and the Couette flow problems are further idealized 
by the assumption that k, and pD are constants. 

The conservation equations for a gas with the above properties can then 
be summarized as follows. ‘The overall continuity equation is unchanged ; 
it is Cp/ct+ V. (pV) = 0. 

Equation (+), the conservation equation for the atoms, becomes 
p[ea,/et+V.Va,]—pDV?x, = b,2,+6, 7 +55. 

The right-hand side, the rate of production of atoms, is temporarily written 

in the general linear form shown. A non-dimensional temperature will be 

introduced and the constants prescribed for each particular problem. 

When 4 is constant the equation of motion is 

p(eV/ct+V.VV) = —Vp+pV?V+in V(V.V). (8) 

In the energy equation (6), the term ¥;h;V, can be simplified to 

> pA: V; = —h, pDV2,—h,pDV2, 


~—(h, —h,)pDVz, 


= Ah*pDV2,, 
when the term [(C,,);—(C,)2]7 is neglected. Hence 
ns as 27 0 2 
p Di = a RV?T + AhpDV7z, +O. (9) 


4. GAS AT REST BETWEEN TWO WALLS 
4.1. Effect of reaction rate 
Before studying the behaviour of the gas described in the preceding 
section in Couette flow, it is instructive to examine its heat transfer 
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characteristics in the absence of motion. Consider, then, the case in which 
the gas is confined between two infinite walls, a distance s apart, at constant 
temperatures 7, and 7°, 

In this situation, the atom conservation equation (4) becomes 

— pD d?x,/dy? = m, K, = b,4,+6, T +s, 
where y is the coordinate normal to the walls. ‘The energy equation (9) 
reduces to 
kd?T/dy? + Ah°pD d?x,/dy? = 0. 

The overall continuity equation and the equations of motion are satisfied 
by V = 0 and p = constant. 

Now write m, K, in the form 

m, Ky = g(%,,9—%), 

where 
T-T,, 
and g and ~,,, are constants. "hus we assume that in chemical equilibrium 
the mixture would be composed completely of molecules at the temperature 
of the cold surface, y = 0, and that at the temperature of the hot boundary, 


é= 


y = 5, the ratio of the density of atoms to the density of the mixture would 
a: %,,- The rate of the reaction is characterized by the magnitude of g. 


To put the equations in non-dimensional form, let 
n= y/S, X= O4/Ains G = (gs*)/(pD), 
_ pDAN%s,, Le Ah x,, 
A(T, 7 Fol al - : 
in which Le is the Lewis number, pDC,,/k. ‘The dimensionless temperature 
6 is defined above. 
In these terms we have 
d*x/dn? —G(«%—6) = 0, (10) 
d?6/dy? + H d*x/dn? = 0. (11) 

The first of these equations is a statement of the balance between the 
atom production in a layer of the gas and the net diffusion out of the layer. 
The constitution of the parameter G reflects the fact that it is the ratio 
of the reaction rate to oe diffusion rate that determines the degree of 
departure from chemical equilibrium. 

In equation (11), H is a measure of the energy carried by mass diffusion 
relative to that carried by thermal conduction. The energy transport by 
diffusion comes about frem the dissociation of molecules in the neighbourhood 
of the hot wall, the diffusion of the resulting atoms towards the cool wall, 
and their subsequent recombination. 

When the walls are not catalytic the boundary conditions on « are 


éx| 


Oo 
cs a = 





_ 


= 0. 





4 =1 


From the definition of 6, 4 )) = 0 and 6(1) 
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Solutions to equations (10) and (11) satisfying these boundary conditions 
and for various values of the reaction rate parameter G are given in figure 3. 
The unbroken line is the atom concentration, and the broken one the 
temperature profile. Only half the channel is shown in this figure since 
the solutions are antisymmetric about the centreline. ‘The parameter H 
has been set equal to 9, an arbitrary value large enough to show plainly the 
effect of dissociation. 





i) 
1 











aig. 
5 


Figure 3. Temperature and atom concentration profiles for non-catalytic walls, H = 9. 
Effect of reaction rate parameter G. ‘Temperature 6 given by broken lines, 
atom concentration « by unbroken lines. 


It can be seen from figure 3 that with increasing G a region in which « is 
almost exactly equal to @, i.e. a region of chemical equilibrium, spreads 
from the centre of the channel toward the walls until finally at large G 
non-equilibrium conditions exist only in thin layers next to the walls. 
Hirschfelder (1957) observes that these thin regions of large change in 
temperature and concentration gradients constitute boundary layers. The 
reason for their existence can be clearly seen in the present simple analysis, 
since the energy transported by thermal conduction and by diffusion outside 
the boundary layers must, when the wall is non-catalytic, be transferred 
to the wall solely by conduction. ‘Thus, the need for a greatly increased 
temperature gradient at the wall is evident. 
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In general, boundary layers will appear in the solution whenever the 
reaction rate is high and the boundary conditions are incompatible with 
chemical equilibrium. 

As G - o the fluid approaches equilibrium throughout the channel. 
‘The non-dimensional temperature gradient, d6/dn, approaches unity except 
at the walls where the limiting value is (1+H). When G - 0 the energy 
carried by diffusion vanishes (since dx/dyj —> 0) and d@/dn > 1 everywhere. 
The concentration of atoms becomes constant across the channel at a value 
in equilibrium with the average temperature, i.e. « = 0-5. 

Figure 4 shows the effect of G on the non-dimensional heat transfer, 
g*, defined by 

gt = aes 

" — R(T,-T,,)' 
where q,, is the heat flow to the wall y = 0 per unit area per unit time. 
This quantity, equal to d@/dn at the walls (when, as in this example, the 
walls are not catalytic), ranges from unity when G = 0 to (1+H)asG => o. 








Figure 4. Dependence of heat transfer on reaction rate, no flow, non-catalytic walls. 


Now to illustrate the influence of a surface catalytic reaction, let the 
reaction rate on both walls be such that the mixture there is held in chemical 
equilibrium. ‘Then the boundary conditions are 

a(0) = 0, a(1) = 1, 

a(0)=0, @(1)=1, 
and the solution is 6 = « = y, regardless of the reaction rate. ‘The non- 
dimensional heat transfer is (1+H), the value at the limit G = o in the 
preceding case. Note that there is a gradient in concentration at the wall 
and hence a flux of energy by diffusion to the wall. 
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Notice, also, that chemical equilibrium exists even for zero reaction 
rate, a result which follows, of course, from the particular boundary 
conditions that were assumed and from the linear reaction rate law. 


4.2. Chemical equilibrium 

This section on the heat transfer through a dissociating gas will be 
concluded by consideration of the simplifications which arise when the 
gas is assumed, at the outset, to be in chemical equilibrium. In this 
discussion the simplifications introduced in §3 are not needed. 

Consider then the energy equation which comes directly from 
equations (3) and (6) for a two-component reacting mixture of gases 
between walls at different temperature 

vd E = tie pi +h, pD | «0, 
dy|_ dy dy 7 dy 
Integrating once and remembering that (z,+,) = 1, we have 


k — + (hy, - h.)pD const. = q,,. 
ay dy 





Now in equilibrium, x, = %,(7,p). ‘Then, since the pressure is constant, 
the above equation can be written 


9D dz, ] dT 
ki 1+(h,—h,) & Lj/—=q, 
| Un) 7 | dy ~ { 


in which the expression in brackets is fixed by the temperature. ‘Thus 


»D dx 
bist ase 
| its) k 7t | 


may be considered to be an effective thermal conductivity, ker. Hirschfelder 
(1957) also notes the possibility of such a definition and remarks that the 
expression is a generalization of the Eucken correction to the coefficient 








the expression 


of heat conduction for polyatomic molecules*. 
Again regarding the mixture as a single gas, we define the effective 


specific heat by 





(C,, Jer = (=) = — (a,h,+ahy) 
/ C 7 » 





‘Then since 


h,= | (C,)dT+Ah?, 


(Cor = (hy —he) = 


- t p’ 


*H. W. Leipmann also made use of this concept in a series of lectures at the 
Douglas Aircraft Company, California, in February and March 1955. 
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where r. equal to [x,(C,,);+%(C,)s], is the average specific heat of the 
mixture when the components are considered to be inert. Now notice 
that if the viscosity is assumed unchanged by the reaction, then an effective 
Prandtl number can be defined by 





| Conf 1 " 8 hg) 4 
(C, erp c. dl 











(Pr)ey = ——Petf — 
; Ret pD dz 
M1+kh-k) 
| la 7? | 
I pe (hy — fy) 4 
_ py cS @ 
1. = hy)Le day . 
Cc, 47 


where the Lewis number is pDC,, k. Hence, when the Lewis number 
is unity (Lees (1956) estimates that for air Le = 1-45), the effective Prandtl 
number of the dissociating gas is the same as that for an inert mixture of 
the same composition. 

Kuo (1957) introduces a constant effective Prandtl number into his 
study of the dissociation of air in the boundary layer. It appears, however, 
that in the transformations of the energy equation which give rise to the 
parameter, the thermal conductivity itself is interpreted as an effective 
conductivity. ‘Thus a direct comparison with the variable (Pr)er above 


is difficult. 


5. COUETTE FLOW 
5.1. Walls at equal temperature 
We turn now to Couette flow, the steady flow between two plane walls. 
One wall is stationary and the other is in uniform motion in its own plane. 
The equations describing such a flow become, for the simplified gas, 
pD d?x,/dy? — g(a, —%, 6) = 0, 
pp. d?u/dy* = 0, 
k d?T/dy? + pDAN® d?x,/dy? + u(du/dy)? = 


where ¢ is a non-dimensional temperature defined below. 


| 


Let us restrict our attention first to flows between equal temperature 
walls and let 





y , u T-T,, o 
”) — . u* =-—-, db = —_—_—_ u P a= pa ; 
s u, jp Can 
Z Dy no 0, 0 
p2 ls Pi Ga § H- pDAhx,, _ Le Ah, 
Cs Fe pD kT, ally il 


where uw, is the velocity of the moving wall, 7’, is the wall temperature, 
s the distance between the walls and Le is, as before, the Lewis number. 
The reaction rate law has again been put in a form which implies that there 


would be no atoms in an equilibrium mixture at the wall temperature. 
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Thus, we have 


° d*x/dyn? —G(«z—¢) = 0, (12) 
d*y*/dy* = 0, (13) 
d*/dy* + H d*x/dn*? + U?(du*/dyn)* = 0. (14) 


‘The solution to equation (13) satisfying the boundary conditions is 
u* =n, 

and the dissipation term in (14) becomes a uniformly distributed heat 
source of strength U*. Hence, for a fixed value of U, the energy flux 
through the walls is independent of G and H and simply equal to the energy 
‘Tt /00_ G=/0 7 
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Figure 5. Couette flow, non-catalytic walls at equal temperature, H = 9. Effect of 
reaction rate parameter G. (a) Temperature profile. (6) Atom concentration 
profile. 

flowing from the source. The effect of the diffusion will then be to reduce 

the temperature gradient (and hence the temperature) required to transfer 

the energy to the walls. 
The boundary conditions are 


4(0) = (1) = 0 


and, when the walls are not catalytic, 


dz. | dx. | ahi 


dy|, 0 as dn , 1 


Solutions to equations (12) and (14) satisfying these conditions with 
(1+H) = 10, U = 1, and for various values of G, are shown in figure 5. 
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In this figure and in the following, ¢,.. denotes the non-dimensional 
temperature at the centre of the channel when the gas does not dissociate ; 
its value is }. The solutions in this case are symmetrical about the channel 
centreline, and again only half the channel is shown. 

Since dx/dyj = 0 on the walls and the energy flow from the distributed 
source is given, the temperature gradient on the walls is independent of G. 
As G increases, again, as in the preceding section, the fluid approaches 
chemical equilibrium except in the neighbourhood of the walls. 

On the other hand, when G —- 0 the temperature distribution becomes 
that for an undissociated gas while « assumes the equilibrium value for the 


average temperature. 
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Figure 6. Couette flow, catalytic walls at equal temperature, H = 9. 
reaction rate parameter G. (a) Temperature profile. (6) Atom concentration profile. 

As in the simple heat transfer analysis, the effect of catalytic reaction at 
the walls will be illustrated by solutions with «(0) = «(1) = 0. Figures 6 (a) 
and 6 (b) summarize these solutions. As expected, a large value of G brings 
the gas to equilibrium. (Note that the « scale in figure 6(b) has been 
changed.) But now as G->0 all the atoms disappear while the temperature 
again becomes that of an undissociated gas. The surface reaction is 
essentially a sink for the atoms. 

The above calculations have dealt with the effects of variation in the net 
Another important property of 


rate of dissociation and recombination. 
This coefficient, D, appears in H 


the mixture is the diffusion coefficient. 





| 
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as well as in G. 
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Therefore, a change in G alone should be visualized as 


a variation in the reaction rate or channel width, not of D. We may vary D 
alone by holding the product GH constant while we vary H. When there 
is no diffusion (H = 0) and the walls are non-catalytic, the temperature is 
unaffected by dissociation (see figure 7) even though the mixture is in 
chemical equilibrium: all the energy is flowing by thermal conduction. 
As D-- w (H -» ~) the atom concentration becomes uniform, as one 


would expect. 
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Figure 7. Couette flow, non-catalytic walls at equal temperature, GH = 900. Effect 


of diffusion coefficient. (a) Temperature profile. (4) Atom concentration 


profile. 


5.2. Comparison with the boundary layer 

\ll the examples just discussed were ones in which the walls were held 
at equal temperatures. We consider next boundary conditions that are 
more suggestive of boundary-layer flows. 1p 
visualized as the edge of the boundary layer and the lower as the wall. 


The upper plate, » = is 

lor instance, the conditions 

A( l ) EC 
dx 


0), 


dn . 0 dy y 0 


would correspond to the flow of an undissociated, cool gas past an insulated, 


non-catalytic wall. 
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The solution in this case is already contained in figure 6 if the centre- 
line there, 7 = 0-5, is taken to be the insulated wall and the station » = 0 
viewed as the ‘free stream’. ‘The fall in the wall ‘recovery temperature ’ 
as G increases is shown by the values at 7 = 0-5. It should be remembered, 
however, that the reason for the reduction in this case is the enhanced 
conduction from the fluid of the dissipated energy. In the boundary layer 
there is an additional effect. ‘The temperature is generally reduced through- 
out the layer because a large amount of the kinetic energy entering the layer 
is absorbed in dissociating the new fluid continually entering. 

Another case of interest is the flow of a dissociated gas past a cool wall. 
The corresponding Couette boundary conditions are 


A(1) = 1, 4(0) = 0, 


a(t )\= Fr, x(0) = 0 or %| = (), 
dy ln =0 
where we have redefined U* and H by 
Bee, Bese, 
Colle fel C,(T,-—T,,) 
and now 
T-T,, 
= ——> 
T,-1 


uw 
where 7, is the temperature at y = s. Equations (12), (13) and (14) are 
unaffected by this change. 
The non-dimensional heat transfer to the cool wall, the quantity of 
greatest interest here, is defined by 
Jw 
Ye <a 
~~ BT ~7) 
and is equal to d@/dy at 7 = 0 when the wall is not catalytic. An approximate 
expression for this quantity, valid when exp{G(1+ H)}!/? > 1, is 
dé | 1 l+H U? U?H 
dy|,-0 1+H{G(1+H)}-? z G(1+A)] 
The heat transfer increases with G and 


lim g* = 1+H+43U?. 








When the cool wall is catalytic, «(0) = 0, and the non-dimensional 
energy flux to the wall by thermal and diffusive transport is 
(d6/dn + H dx/dn), — 
a quantity equal to 1+ H+ 4U? whatever the value of G. 
If the gas is assumed to be in equilibrium initially, then « = 6, and 
equation (14) (with ¢ replaced by @) yields 
6+ HO = —}U*%y?+C,7+C,. 
The conditions 6(0) = 0, @(1) = 1 require 
C, = 0, C, = 1+H+3V?, 
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and thus 


Pee ee =1+H+iU". 
dy |, 0 


wn 
as 


A mechanical analogy 
A mechanical analogy, helpful in visualizing the form of the solutions 
discussed in this and in the preceding section, can be developed as follows. 

Consider two strings, one weightless and under tension H, the other 
weighing UU per unit length and under unit tension. Denote the weightless 
string displacement by « and the other by 6. Let the strings be connected 
by a distributed spring, an elastic sheet, with spring constant per unit 
width equal to GH. 

The force balance on the « string requires 

H d*a/dyn? = —GH(@—-x), 
and the @ displacement satisfies 
d*6/dy? = — U?+GH(6—-«). 
A rearrangement yields the Couette flow equations and setting U? = 


QC 


vields the simple heat transfer case. 





A | 
, f 
1 j 
; L/] f 
7 : : 
; a @ 


j 2 
U 


Figure 8. Mechanical analogy. 











In the sketch (figure 8) the boundary conditions are those for catalytic 
walls at equal temperature. We can enforce the non-catalytic condition, 
(dx/dy)wai = 0, by connecting the « string to frictionless sliders on the 
wall. ‘lhe vertical component of the force applied to a wall by the strings 
are dé/dy and Hdxz/dy. ‘Thus, these forces are analogous to the energy 
flowing into the wall by thermal conduction and diffusion, respectively. 
A variation in the reaction rate constant g corresponds to a change in the 
spring constant GH with all other parameters (including H) fixed. 


6. LAMINAR BOUNDARY LAYER 


6.1. Solution of the energy equation 


In this final section we discuss the flow of a gas past a flat plate in 
circumstances in which dissociation and/or recombination may occur 
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within the boundary layer. ‘Three recent papers have dealt with the heat 
transfer in laminar boundary-layer flows of a dissociated gas. Lees (1956) 
studied two limiting cases of heat transfer to blunt-nosed bodies at 
hypersonic speeds, the case of chemical equilibrium, and that of no 
recombination except possibly at the surface within the boundary layer. 
Fay & Riddell (1956) obtained affine solutions in the neighbourhood of 
the stagnation point for those limiting cases and for a finite recombination 
rate. Kuo (1957) studied the laminar boundary layer of a gas in chemical 
equilibrium on a flat plate. 

We again consider the gas to consist of two components, atoms and 
molecules. ‘The usual boundary-layer approximations* to equations (2), 
(4), (5) and (6) in §2 are: 


< (pu) + = (pv) = 0, (15) 
C oy 

o( ee z) - = (oD) == m, K,, (16) 
Ox Cy oy 


(us, to oe) = =[% oF (hy —hy)pD =| +o). (18) 
Ox oy oy oy oy oy 


/ 


Here the pressure has been assumed constant, which is, however, a poor 
approximation near the leading edge of a flat plate at velocities high enough 
to dissociate air. The equations might be better associated with the flow 
past a surface shaped to minimize the pressure variation, a wind-tunnel 
wall, say. In the analysis which follows immediately below, the simplifi- 


4 


cations introduced in §3 are not necessary, and they will be first introduced 
in § 6.2. 


Lees (1956) has shown that the energy equation takes a much simpler 
form when the Lewis number is unity. ‘To see this, note that 


h = ah, +h, 
with 
Then oh = OT Ox, 
— p Aa. 


oy ( 





where C,, = 4,(C,,)1 4 x(C’,)o. Thus, the first term on the right-hand side 


* Reaction rates may in many cases be high enough to confine the region of 
departure from chemical equilibrium to the neighbourhood of the plate leading edge. 
The use of boundary layer equations to calculate the non-equilibrium flow under 
these conditions would, of course, be incorrect. 


12 
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of equation (18) can be written 


=| = oh + (h, —hy) & (Le—1) a 
cy c: oy Be cy 


where, now, Le = pD( a k. Hence, when Le = 1, the term 


contains the energy flux due to diffusion as well as the usual heat flow 
caused by the thermal gradient. Now the energy equation can be written 


o( u = +Y =) = — (= =) + nS) : (19) 
CX oy oy fe oyV oy 


If we have boundary conditions on / or conditions on 7 and ¢, that determine 
the boundary conditions on A, then this equation together with (15) and (17) 
may be solved by the procedure of Chapman & Rubesin (1949). ‘This 
solution is obtained without reference to equation (16) and, hence, does 
not require knowledge of the reaction rate. ‘The individual distributions 
of temperature and atom concentration, however, would not be known. 
Interesting flows in which the enthalpy boundary conditions are known 
at the wall are flow past an insulated and non-catalytic plate, for which 


in) 
— 


1 — oT Je 
= ( p a T (h, tas hy) bie = 0, on y = 0: 
oy cy cy 


XQ 


and flow past a wall at known temperature on which a catalytic reaction 
maintains chemical equilibrium, for which 
h(x, 0) = «,(7,,)h,(T,,) + %0(T,,)ho(T,,), 

in which the «’s depend only on temperature since the pressure is constant. 

We assume that the properties at the edge of the boundary layer are 
known and are independent of x. Again this is not a realistic assumption 
for the flow immediately following a strong leading-edge shock. ‘The 
Chapman & Rubesin procedure then follows. First introduce the stream 
function defined by 


pu =p Cus ey, pv = — px Cus Ox. 


Then in x, % coordinates the momentum equation becomes 


Cu 1 a Ou 
er See aoa b il u po . 
ox Px ¢ ub Ps ous 


Next we assume pp = Cu,p,, with C a constant, then 








ou Ch, @ _ Cu 
Ox * Px Cus ( =) 

Lees (1956) and Kuo (1957) justify this approximation for air in chemical 
equilibrium at high enthalpy levels. Here, where equilibrium does not 
necessarily exist, the approximation, even if less accurate, has little influence 
on the nature of the following qualitative results. ‘The energy equation (19) 
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becomes, when the Prandtl number C,, / is constant, 


oh ypxC Od ‘ oh Ce. . eu\? 
ox p. Pr os Cus Px ne , 











Let 
oe oe pe = Poke ae ae 
u, h, fs 
L Ps [(Mxc/px ux LC]'” 


where L is a reference length and the subscript oo denotes free stream 
conditions. ‘Then we have 





, ‘ ae 
De siecle (w a (20) 
ox* ous* Ous* 

and 

om  t @ , oh* ux. 5 pe 2 (21 
ax* ap" at ag =) on) 


The solution of (20) is 
u* = if'(n), 
where 7 is defined by f(y) = %*/x* and fis the Blasius function satisfying 
if +f" =0, 
f'(c) =2, f'(0) = f(0) = 0. 
Since u* is a function only of 7, it is convenient to write the energy equation 
in x*, » coordinates 





ozh* . ch* , , oh® Ba 
. 3 + Pr f — —2Prf'x* seal = —1}Pr hi 7’ ¥. 
on* wi | : ox* h. : 


Chapman & Rubesin find, by separation of variables, a solution of this 
equation (in their analysis, 7 is the dependent variable) for a wall 
temperature that varies with x. The solution will, of course, be applicable 
here when the wall enthalpy is a given function of x. 

When the wall is non-catalytic and insulated, the boundary conditions 
are h*(x*, ©) =0 and ch*(x*,0)/e7 = 0. These conditions are satisfied 
by a solution that depends only on 7». ‘This particular solution, denoted 
by N(), satisfies 

N” + PrfN’ = —4Pr(u2,/h,)(f’), 


from which we get 
N = 3(u2,/h,)o(n), 


where ; 
o(n) = Pr | LPP | Lf Wy P-P doa. 
Jn /~ 0 
The function o(7) for Pr = 0-72 is given by Chapman & Rubesin (1949, 
figure 4). 
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If the free stream is undissociated and the gas in the boundary laver 
assumed to be in chemical equilibrium, then the above solution becomes 
essentially that obtained by Kuo for this case. 


6.2. The temperature and concentration profiles 

As remarked earlier the above solution yields no information about the 
temperature and concentration profiles when the reaction rate is finite. 
The qualitative behaviour of these variables can be obtained if we introduce 
approximations of the kind discussed in §3. Again let h = C,, 7+ Ah°x,, 
and assume that 

m, K, = (p/p. )g(%,, 7* —%), 

where 7* = (7 —T,)/T,. ‘Thus we are assuming that there are no atoms 
in the free stream and that any increase in temperature above the free stream 
value causes dissociation. ‘lhe factor p/p, is introduced into the reaction 
rate law for convenience; when it is present the energy and concentration 
equations are linear in the x*, 7 coordinates defined above. 

With these expressions for h and m, A, and with the help of equation (16), 
the energy equation becomes, in the x*, coordinates, 








a2T* .ol* oi ion™ Pr u> aa 
——— Pr f —— - 2Pr f y* —_ ——— (f )-+ 
ene cn ox" 4¢ p fr. ? 
4Pr LAh°g(a,,, T* — x,)x* 
patna ts 
and the concentration equation is 
at led | On mye OO 4ScL 2(a,, T* —a,)x* 
1 4 Sof 1 -2Scf'x* 3 = AG =. 
( | ( ” ’ ( x* Ps u, 


where Sc = /pD is the Schmidt number. We see from these equations 
that the significant x variable for changes in 7* and ~, is 


* Log vg 
pu, p,Uu, 


Calling this quantity x’ and letting x = ,/z,,, we get 








a2T* _oT* , ,ol* 
—— Pri pene 2 2Prf x9 = 
( n~ P ( ” . Cie 
Pris 8 4Pr Mh°z,,, et _ x) (23) 
4C,T, * i? 
—— +Sef — —2Scf'x®— = —4S8ex"(T*—«). (24) 
cn” coy) ON 


Now, following Marble and Adamson (1954)* we could write 
T*( x. ») < T*(n) 1 TH(y = 4 T* (ny dae" 


a(x®, ) = 4,(7) +4 X11 (9 )x° t. o111(9) x 4 ae 


* Equations similar to (23) and (24) arise in a flame propagation problem treated 
by these authors. Howarth (1938) formulated a series solution of the form shown 
for an incompressible boundary layer problem. 
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If these series are substituted in equations (23) and (24) and like powers 
of x” are equated, we get 
— 
Peas 


4C,, I. 


p 


T}" + Prf Ty’ = - (f") 


4Pr Ah’x,, 








T* + PrfT*' —2nPrf'T* = a (T* -%,-1); 
‘ me as 3% , 
and 
oy : Sc fx, = 0 
a" + Sc fa, —2n Scf'x,, 4Sc(T*_, — %n_1)- 


Thus, the coefficients could be determined term by term when the 
boundary conditions depend only on 7. Numerical or machine solutions 
of the successive equations for 7* and «, would be required, however, and 
in view of the approximation to the reaction rate, are not justified. 

We note, however, that for the case of an undissociated gas flowing 
past an insulated and non-catalytic plate, i.e. when the boundary conditions 


are 
aT *(x° (< 
T*(x°,0)=0, =O) 6, 
ov) 
a 0 ( 
x(x, 0) = 0, da(x?, 0) ih, 
oy 


then the leading terms are those that would be obtained for an undissociated 
gas; 4, is zero and 7* is the usual temperature distribution in a perfect gas 
on an insulated plate. A heuristic argument, given below, suggests that far 
downstream the temperature profile returns to this shape (but with a lower 
absolute magnitude). We may expect from this that the change in the 
temperature profile shape is relatively small and that the use of similar 
profiles in integrals of the energy and concentration equations would yield 
approximate equations for the change in magnitude of the temperature and 
concentration with x. 

To see why the temperature may be expected to have the same shape 
far downstream as at the start, recall that in Couette flow an increase in the 
channel width increased the reaction rate constant and brought the gas 
closer to chemical equilibrium. Likewise, here the concentration gradients 
across the boundary layer which tend to keep the gas out of equilibrium 
become weaker as the thickness increases. ‘Therefore, since gradients 
along streamlines also decrease with x, it is reasonable to expect that the gas 
approaches chemical equilibrium as x increases. We can easily find the 
temperature and atom concentration at equilibrium if we assume that the 
Lewis number is unity. First, we have the solution for h* in this case; 
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it is the function N(7) obtained earlier. ‘Then, since we are assuming that 
h = C, T+Ah®a,,% and (the equilibrium condition) « = T*, we have 


( 
en ee (14 o)t 





h, Cor. 


p 
Therefore, in equilibrium, 


T* =Non/(1+ 3 " Ah' a3) 


Now observing that N and 77 satisfy i same equation and boundary 
conditions, we see that the temperature profile has, in equilibrium, the 
same shape as it has at the leading edge but, of course, is reduced in 
magnitude. ‘Thus, if, as the above considerations suggest, the gas approaches 





equilibrium in the boundary layer as it moves downstream, we can describe 
approximately the approach to equilibrium by setting 
T* = N(n)Q(x°), — & = N(n)A(x°), 
in integrals of the energy and diffusion equations (23) and (24) to determine 
Q and A. 
We have then 
Q | N’dn+PrQ| fN'dn—2Prx°Q'|  f'Ndy 
~ 0 ~ 0 *~ 0 
PFU: (CO ome 4Pr Mh®x,,, x° *o 
—~ gee | (Pent — ae 
Totats hy 
and 
A| N’dn+ScA| fN'dy—2Scxd'| f'Ndy 


~ 0 0 “0 


-4Sex(Q-A)} Nady 


Let 
b= gig [Uh 
and note that since 
| N" dy = 

equation (22) shows that : 

=-—| {Nd 
Further ™ 

| fN'dy=—[ f'Ndy=—-], 

Then, if we let 

I,=4[ Ndn 
we have finally : 

— Aha, Ty 


¥9(Q2—A) (25) 


‘a: a 3 
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and 


A+2x° A’ = > x°(Q—A). (26) 
1 
Eliminating Q and setting 


0 
$= a(1 A - HA, 





7, 
m 1 1+ AA°x,,, I, x? 
ile 2 C, yl 1 
we get 
276" + 2(1 +2)8' + 3(z—3)d =s, (<7) 


where the derivatives are taken with respect to z. The initial conditions 
are Q(0) = 1 and A(0)=0. (It is interesting that any other conditions 
at x° = 0 in equations (25) and (26) result in infinite derivatives in Q and A 
there.) The condition on 2 can be replaced by a second condition on A as 
follows. Differentiation of equation (26) yields 


A’ +2x°A" +.2A’ = (I,/1,)[x9(Q' — A’) + Q—A], 


from which we get, by the provisional assumption that derivatives at x° = 0 
are finite, 

: dA 1 

A’(0) = — (0) =- 

dx® 3 


I, 
| 
In terms of 6 then, the initial conditions are 

d3(0) _ 


az 


6(0) = 0, 6'(0) = 
Solutions to the homogeneous equation corresponding to (27) may be 
found in Kamke (1948, p. 451). They are z~? and z-!?e-*. From these, 
a particular solution to equation (27) satisfying the required boundary 
conditions may be found by a formal application of standard methods. 
The solution is 

6=2-2 va tel 1/2e! dt. (28) 

~~ 9 
An expansion of this expression about z = 0, or, better, a series solution 

of the original differential equation, yields 


@ 
o= > ear. 
i 0 
with % = 4/3 and 
a, 
ea x en) eee 





In figure 9, the solid curve extending from z = 0 to s = 4 is a plot of the 
first few terms of this series. 

The behaviour of the solution at large z also may be determined either 
from an examination of the integral in equation (28) or from a series 
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formulation. Erdeélyi (1956, p. 69), states that the integral 


| ett—” dt 


Jo 
behaves like 
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Figure 9. Atom concentration level as a function of distance along the plate. 


The use of this asymptotic expansion in the expression for 5 in 


equation (28), or the substitution of a series in inverse powers of z directly 
into the differential equation yields 


6= VY¥bx2- 
— n ’ 


U 


where 5, = 2 and 


b, .1 = (n—})b,, n=0,1,2,.... 
Thus, as 2 


~ ¢ 


ipproaches infinity, 6 approaches 2, and since 


A = 48(14 Amn)" 
i 7. 


p 
we have 


lim A - (1 , Aen 
oF. 


p 
Equation (26) shows that 2 approaches this same value. 


Remembering 
that 


4. T* = N(n)Q(x°), 


x = N(n)A(x°), 











Dissociation in Couette and boundary-layer flows 139 


we see that the fluid approaches a state in which 7* = a, i.e. chemical 
equilibrium. 

The above series is divergent. ‘The values of 5 shown in the dashed 
curve in figure 9 were computed by retaining, at each value of z, the term 
up to but not including the smallest. (See Jeffreys & Jeffreys (1950) for a 
discussion of this procedure. ) 

The concentration of atoms, then, rises from zero to 

Non(t + ct) 
ay 
while the non-dimensional temperature falls from N(y) to this same 
quantity. 





The author wishes to thank Professor H. W. Liepmann for his 
stimulating and encouraging interest in the problem. He is also grateful 
to Dr Gilles M. Corcos who read the original draft and suggested many 
improvements. 
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SUMMARY 

‘The approximate T’ricomi equation relevant to sonic speed of 
two-dimensional small-disturbance flow is solved by separation of 
variables, where these are certain stipulated functions of the 
Cartesian velocity perturbations. ‘The symmetric flow patterns 
obtained from this solution are shown to correspond to those 
about half-bodies, whose ordinates vary as x” where 0-4<n<1, 
and x is the distance along the plane of symmetry. ‘The surface 
pressures on such bodies are deduced. 

In particular, the body whose ordinates vary as x? has a sonic 
surface velocity, except at the nose, where an edge-force which 
causes a drag force is shown to exist. On the assumption that a 
free-stream breakaway (at sonic velocity) occurs at the shoulder 
of a body, this solution thereby yields the flow about an aerofoil 
of the same shape having a flat base. ‘This bluff-nosed section 
has only a little more than half the drag of a wedge on the same 
chord and base. 


INTRODUCTION 


The use of the hodograph transformation of the equations of inviscid, 
transonic fluid flow is too well known to need description here, particularly 
in its approximate form which leads to the Tricomi equation. Whilst 
this is not the place for a rigorous criticism of the assumptions involved 
in this approximation, it will be recalled that its essential basis is that the 
fluid perturbation velocity is in general of small magnitude compared with 
the speed of sound. It is not unusual for this approximation to be used 
fora flow about a body witha stagnation point, though plainly the assumptions 
made break down in the immediate locality of such a point. Indeed the 
perturbation velocity becomes singular at this point. However, by reducing 
some relative proportion representing the body geometry and governing 
the magnitude of the disturbance created—such as the (thickness/chord) 
ratio of a wing section—it is possible to arrange that the perturbation of 
velocity exceeds any prescribed limit only within a region whose dimensions 
are small compared with those of the body. 

In the present paper we shall discuss the flow about b/uff bodies, for which 


there is inevitably a ‘stagnation point’ represented by a singularity in the 
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perturbation velocity, whose form differs according to the shape of the 
body. For instance, the rise in speed along the body surface from the nose 
will be shown to be more rapid on bluff shapes than on those with more 
nearly sharp noses. However, we are also dealing with half-bodies, whose 
scale is characterized by a single length—such as, for example, the nose 
radius of the body if it happens to be parabolic. It transpires that the 
approximation is inapplicable within a region about the nose whose 
dimensions are of the order of this characteristic length. ‘Thus the solution 
is only valid as applied to the flow at such distances from the nose (large 
compared with this length) that the perturbation velocities are bounded 
within prescribed small limits. In particular, the surface conditions are 
only correctly simulated where the surface slope is small. Relative to this 
scale, the region of failure is small, and by the same token the actual shape 
of the body surface in this region is really immaterial. 

A precise distinction is drawn if it is stated that, by using the Tricomi 
equation, we can strictly deal with only the asymptotic behaviour (towards 
infinity) of the sonic flow about certain half-bodies, whose asymptotic shape 
is stipulated. Asymptotic conditions are approached within any desired 
accuracy at distances which are large compared with the characteristic 
length of the body asymptote (which we shall later identify by the symbol e). 
Usually a solution is sought of this Tricomi equation by the method of 
separation of variables, these being in fact the Cartesian perturbation 
velocity components of the fluid motion. We seek here a solution of the 
equation by the same method, using, however, two different independent 
variables which are functions of both these velocity components. These 
particular functions have appeared before in various essays on the subject 
of transonic flow, though it does not appear that they have been used in 
the way employed here. 

It may be that, by superposition of such solutions as those obtained 
below, it is possible to determine the flow about certain closed bodies; 
it is certain that various asymmetric flow fields about lifting surfaces may 
be derived in this way. ‘These could be of practical interest, though the 
appearance of supersonic regions of flow and limit lines greatly complicates 
such analysis. However, we shall content ourselves here with the simple 
basic solutions which yield the flow about semi-infinite and symmetric 
bodies whose ordinates vary as x”, where x is the distance downstream of 
the nose in the free-stream direction. Provided the power index n is between 
0-4 and 1, the flow velocity is nowhere supersonic. The half-body whose 
ordinates vary as x” is of particular interest since, within the accuracy of 
the approximation, the surface speed is constant and sonic; because of 
this the flow may also be envisaged as that about a finite aerofoil of this nose 
shape having a flat base, with free-streamlines extending from the corners 
of its base to enclose the wake. It might be concluded that such-an aerofoil 
has zero pressure drag, but this is not so because of the existence of an edge 
force at the nose; however, its drag is much less than that of a wedge on 
the same base and of the same chord. The parabola is another member of 
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the family of half-bodies which is also of interest, since it coincides with the 
displacement surface of a laminar boundary layer on a flat plate. The 
semi-infinite wedge does not yield a bounded field of velocity and is not 
amenable to treatment. 

In general, we treat only the surface conditions and not the complete 
streamline pattern, though computation of this would be possible. ‘The 
hodograph solution is first obtained and then its interpretation in the 
physical plane is discussed. Finally, the particular applications already 
mentioned are amplified. 

If a liberal view is taken, it can be suggested that the results lend support 
to the general belief that a bluff-nose can provide less drag than a sharp 
one at the speed of sound. The results may prove of use in design problems 
of both internal and external aerodynamics; and they enable the extent 
of the inviscid wake to be calculated in conditions of shock-free flow 
separation at sonic speeds. 


‘THE HODOGRAPH SOLUTION 
Suppose that the velocity of flow of a gas is everywhere nearly sonic; 
then the Cartesian components U and V (with the direction of the free 
stream in the direction of the positive x-axis) can be written as 


U=a(i-—2), are wk, (1) 


where wu and @ are small compared with unity, and a* is the critical speed 
of sound. Using u and wv as independent variables, it can then be shown 
that the stream function ‘’ and velocity potential ® of the flow obey the 








equations 
®,= u¥,, ®,=-¥ 


’ v ue 


‘The second of these relations is satisfied by defining a function (Q(u, v) such 
ii a 

that ®= -Q,, P=Q,, (2) 
and the first then yields the well-known Tricomi equation 

C2 +uQ,, = 0. (3) 

Solutions of this equation by means of separation of the variables u 

and v have frequently been found of interest. We here employ the same 

means, but we separate instead newly introduced independent variables, 


w= 1—+(u3/w). (4) 


We shall concern ourselves only with flows without supersonic regions so 


2 Se on 


w = =u°+ Vv", > = 


that u>0; hence w>0, and 0O< €< 1. If supersonic regions were 
present the value of € would not lie in this range. Positive and negative 
values of v are distinguished by the sign of €'*._ Using these new variables, 
equation (3) becomes 


wt 1wQ,.+ (} -— 7€)Q.+ €(1-€)Q.- = 0. (5) 


uw 6 


If we seek a solution of the form 


Q = f(w)g(E), (6) 
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then the functions f and g satisfy the equations 
wif + Jef’ + GH = 0, 
é(1—€)g” + (4—2é)e’ —(—-n2)g = 0, 


where, for convenience, (4,—") has been chosen as the constant of 


) 


separation. 
The first of these equations may be solved to give 
f(w) = (C, aw" + Cw)! 2, (8) 
where C, and C, are disposable constants. ‘The second may be reduced 
to Legendre’s equation, although its solution in the range € < 1 may be 
more conveniently expressed formally in terms of hypergeometric functions, 
g(é) = BF +0 b—w3 33 2) + Bee FR E+u gw 4:8, (9) 
where again B,, B, are adjustable constants. 
Since the hypergeometric functions are not affected by the sign of p, 
it is no loss of generality to put, say, C, = 0 in (8), and take the value of C, 
as unity. Further, we shall restrict the analysis to flows symmetric about 
the free-stream direction so that ‘is an odd function, and ® and Q are even 
functions, of v (i.e. of €'?). Thus we take B,=0 and it follows on 
substituting from (8) and (9) in (6) that a typical solution of the form 
required is 


_ gn#—112 FY 1 1 e «-€ 
O= Bie F(i+h,5—13 4; €). (10) 
Using (2) and (4), together with (10), we obtain after some rearrangement 
1 > 9 ” 
—D= — 3u’w-™ F(m,£—m; 1; &) 
a*l S ao 
= — (3)!3qy-m+23(1 — £)23 F(m,i—m; 1; €), 
: | (11) 
— Y = 2(m— 3)vw-" F(m,?—m; 3; €) 
a*| ' : 
= 2(m— 3)w-™ t2Eu2 F(m, + —m; 3; &), 
a ne < r7enience re have sas ae : anlacre 1 
where, for convenience, we have put m = 3—p, and replaced 2(u—4)B, 


by a*/ > 0, say. 


THE INTERPRETATION OF THE FLOW IN THE PHYSICAL PLANE 

To interpret the flow pattern in the physical plane determined by the 
solution (11), we need to invoke the following known property of the 
hypergeometric functions, which are of course equal to unity at € = 0 
and are bounded and continuous in the range 0 < € < 1, except possibly 
at €=1. Provided that m> 3, there exist numbers €,, € where 
0<& <& <1, such that F(m,i—m;3;€)>0 for 0<&<&, and 
vanishes at €=&,; and such that F(m,:—m;};£&) changes sign just 
once in 0< & < &, at €=€,. If m< 3 it may be remarked, in passing, 
that supersonic regions of the flow and limit lines make their appearance. 

Interpreting the powers of w by their positive real values (and supposing 
m > 3) we see that at infinity (where ® and ¥F are infinite) we have w = 0 
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(and € arbitrary), so that from (1) and (4) the velocity has there decayed 
to that of the free stream. Further, the point ® = ‘¥ = 0 corresponds to 
a singularity in velocity (w= x —_ € arbitrary). We see that ‘’, but 
not ®, changes sign with v (and €'*) as required earlier to meet the 
provisions of flow symmetry, so that we may restrict ourselves to consider 
the upper half-plane where ‘¥ > 0. For all negative ®, this region is bounded 
by the line ‘* = 0 on which € = 0 but w takes all (non-negative) values, 
that is, on which wu is positive but v = 0. Evidently this is a ‘stagnation’ 
streamline, the ‘stagnation’ point at ® = ‘Y = 0 being a singular point 
of the velocity distribution, as is the usual consequence of the approximations 
involved. ‘The line € = €, then maps point by point on to the equipotential 
line ® = 0 extending from the singularity, the region 0 < € < &€, for all 
non-negative w transforming in a one-to-one relationship to all points 
within the region upstream of this line. Finally, we note that the line 
€ = &, defines the continuation of the stagnation streamline Y= 0 with ® 
taking all positive values, the region €, < € < €, mapping on to that down- 
stream of @ = (0. Since €, and so also v, # 0 on ‘Y = 0 downstream of the 
stagnation point, the stagnation streamline divides at ® = 0 into two 
branches (symmetric about the upstream stagnation streamline), the value 
of v being equal and opposite on each side. We see that the flow is a 
symmetric one about a half-body extending downstream of ‘Y = ® = 
with no singularities in the flow external to it. 

The shape of the half-body in the physical plane will now be determined. 
In so doing, we will suppose that the value of €, = &,(m) has been found 
(as already defined) to be the least positive real root of 


F(m,=—m, 3; &) = 9, for m2> 3; (12) 

for brevity it is convenient to put 
(3)'3(1 — &,)?3 F(m,i—m; 1; €&) = —6/l, say, (13) 
where 6 replaces / as an arbitrary positive length, (6//) > 0, and is bounded 


for any m > 3. We see then from (4) and (11) that on the upper surface 


of the half-body (‘F = 0, ® > G, v > 0) 
v = (€,@)*?, > (14) 


all powers having the positive real interpretation. ‘To identify the velocity 
field with the Cartesian coordinates (x,v) of the physical plane, we note 
that on a streamline changes of position coordinate are related to changes 
in ®, within the accuracy of the approximation used in deriving the basic 
equations (2) and (3), by the expressions 


dx = d(® a’), dy — _ d(D a*). (1 


‘ma 
— 


If the stagnation point ® = '¥ = 0) is the origin of the Cartesian axes, it 
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follows on substituting from (14) in (15) that on the upper body surface 
(f = 0) 





whence on integration, its shape is given by 


(y/e) = (x/e)", 


6m—7 eo" G—n) } (16) 
where n = ——— and «=| ——*— é. 
6m — + n(y +1) 


Since m > 3 it follows that the index m can take any value between 0-4 
and unity. 
The pressure coefficient C, is related to the velocity changes with 


sufficient accuracy by the equation 





s 
c. re (17) 
whence from (14), (15) and (16) it follows that on the body surface 
Cp = 2(n)(dy/dx)?® = n?8 a(n)(e/x)?0-8, 





8 cal 1/3 (18) 
where a(n) = | (- r =) | ; 


The function x(”) may be determined by solving (12); this demands 
a knowledge of the behaviour of the quoted hypergeometric function, which 
may, however, be related to evaluated transcendental functions if the value 
of m is an integral multiple of 1. In particular, we note that for m— 3+0 


2 5 8 \13 x | 
n>—=+0, &—-1-0, x(a) ~ 7 ( . B a (5n— 2), 
) : 1I8\y+1 a2 | 


and form>+ca (19) 


P 7 72m? 13 
n -r— oe. fo ~ — , x(n) ~ neogerarsaee ; 
~  4m* m(y +1) 


where B(a, 4) is the beta function. 

Thus » varies from 0 to ~ as m varies from 0-4 to 1. Its variation is 
sketched in figure 1, composed from a discrete number of solutions of (12). 
The expressions for the hypergeometric function appearing in (12) with 
values of 6m between 4 and 9 are summarized for convenience in table 1. 
The expressions with larger integral values of 6m both for this and also 
for the other hypergeometric function of the fundamental solution (11), 
may be determined from the tabulated functions by elementary operations. 
The manner of these operations is indicated in table 1. 

The point made ‘n the Introduction about the breakdown of the theory 
close of the stagnation point is best appreciated by referring to equation (18), 
from which it will be seen that C,, and dy/dx are in general finite on the surface 
for finite (x/e), where (from equation (16)) it will be observed that (y/e) 
is also finite. In general the velocity perturbations will only be small, as 
required by the Tricomi approximation, at distances from the origin which 
are large compared with e. Thus if it is required that Cp (or u) has an 
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Table 1. Closed expressions for the hypergeometric functions F(m, {j—m; 3; €). The 
*s I = ™ 


expression Be(a, 5) | 1¢-1(1 —t)b -1 dt is the incomplete beta-function. 


order of magnitude (1/.\)), then v has a smaller order of magnitude (1/.V*?), 
and it may be shown that this applies in general where 


, é 


The singularity representing the stagnation point becomes more intense 
as the degree of bluffness, represented by the exponent (1 —v), increases. 
Thus equation (18) shows that as n is reduced towards 0-4 the surface 
speed recovers more rapidly from the ‘stagnation point’, until in the limit, 
on the body shape given by n = 0-4 it jumps discontinuously from infinity 


Some typical examples are illustrated by figure 2. 


to the sonic speed. 
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Figure 2. Comparison of nose-shapes on some of the half-bodies, together with the 
corresponding pressure distributions and drag coefficients on the frontal area 
(over portions shown). Scaling the ordinates by a factor (1/N) reduces the 


pressure and drag coefficients in the ratio (1/N)?/8, 


SOME PARTICULAR SOLUTIONS OF INTEREST 
The simplest case n = 0-4 (or m = 3) is of particular interest, since it 
will be seen from (16) and (18) that it yields the flow past a half-body 
(whose ordinates vary as x”°) over the whole of whose surface the velocity 
is sonic and the pressure equal to that of the free-stream. It might therefore 
be thought that the drag on the body is zero, but an ‘edge force’ exists at 
the stagnation point, as may in fact be demonstrated by considering the 
drag derived from the pressure distribution along any complete streamline 
other than ‘’ = 0. It is more simple to show this, however, by noting that 
the drag on that part of any of the bodies defined by (16) ahead of the plane 
x = const., where the pressure coefficient is C,, can be expressed by means 

of a coefficient C,, based on the local frontal area and given by 

3 
Cy = 5 Co. (20) 
Thus in particular, allowing n + 0-4 from above, we find from (20) with the 
help of (16) and (19), that the drag of the half-body vy = «*°x*> ahead of 
the plane x = const. > 0) is 

5(5)'°B(h, B)spta*e(y + 1). i 
This is independent of x and so the drag originates from the nose. The 
type of singularity in the hodograph which yields such an edge force is 


found by placing m = 3 in (11). 
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This appearance of an edge force is interesting, and it seems to be the 
first occasion that it has been noticed in the study of transonic flow. A finite 
force associated with a singularity is known to appear at a radiused leading 
edge in the treatment of the flow about it by linearized theory, provided 
only that the free-stream velocity component normal to the edge is subsonic. 
\lthough rounded leading edges are strictly outside the scope of linearized 
theory, it is further known that nevertheless the forces are correctly given. 
It is not possible to draw the same conclusion in the present instance, 
because there is no more exact treatment available of particular cases by which 
the issue may be judged. ‘The most that can be said with assurance is that a 
body whose shape is asymptotic to a curve y = e? °x?° has a drag only by virtue 
of the pressure distribution overthesurfaceat distances of order e fromthe nose. 

If the hypothesis is adopted that a free-streamline extends from the 
shoulder of a flat-based aerofoil, enclosing fluid at the free-stream pressure 
(so that there is no base-drag), we see further that this same solution (with 
n = 0-4) then yields the sonic flow about a finite aerofoil (whose ordinates 
y « x25) with the free-streamline forming a continuation of the surface 
downstream of its flat base. ‘The ultimate width of the ‘wake’ contained 
within the free-streamlines would be infinite (at infinity downstream), 
which is the same conclusion that is reached in examining the flow about 
a flat-based wedge. ‘The drag coefficient (on frontal area) of the curved 
aerofoil thus deduced is, from (16) and (21), 

Cy say 1:20( Vmax cy 3(y +i 
where c is the chord and ymax the semi-ordinate at the base. For a wedge, 
the numerical factor is 1-89, instead of 1-20 (Helliweli & Mackie 1957). 

If this same hypothesis about the wake is made in dealing with flow 
about an arbitrary aerofoil, with separation at the point of maximum 
thickness, we see that the shape of the free-streamline representing the 
wake boundary will have ordinates asymptotic to y = e°x?°, Further, 
from a consideration of the fux of momentum at infinity, it will be evident 
that the characteristic length € of this asymptote is determined by the drag 
of the aerofoil in accordance with the expression (21). These deductions 
are in fact confirmed by the work on the wedge previously mentioned. 

The half-body given by y « x!? is also of some interest as, on account 
of its parabolic shape, it evidently shows the displacement effect at sonic 
speed of a laminar boundary layer on a semi-infinite flat plate aligned with 
the stream direction. ‘The self-induced pressure field of the boundary 
layer is then seen from (18) to vary as 1/x'3. The theory breaks down 
when applied to the semi-infinite wedge (m= 1) since it predicts an 
unbounded surface velocity. ‘This is perfectly compatible with results 
for the finite wedge, because these show that the surface velocity becomes 
infinite at the nose; the semi-infinite wedge, on the other hand, has no 
characteristic length, and so the surface speed must be constant, with a 
value corresponding to that at the nose of the finite wedge. 
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SUMMARY 


This paper describes an experimental investigation of the form 
of the large scale motions in turbulent flow. ‘These motions have 
been found to be more ordered than has usually been supposed and 
their origin and dynamics are discussed in terms of physical models 
of typical eddies. 

Nine components of the double velocity correlation tensor 
have been measured at a number of positions in the wake of a 
circular cylinder and in a ‘flat plate’ boundary layer. ‘These 
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have been supplemented by measurements of correlations with 
separations in directions other than the axial ones. In the wake, 
the correlations at large values of the separation are explained in 
terms of two types of large scale motion. One of these is a pair 
of vortices, side by side and rotating in opposite directions with 
axes aligned approximately normal to the plane of the wake. ‘The 
other typical motion is a series of jets in which turbulent fluid is 
projected outward from the core of the wake. It is suggested 
that these are the result of an instability of the turbulent shear 
stress. A qualitative explanation of the apparent structural 
equilibrium of the wake is given in terms of this instability. The 
vortex pair eddies were not found in the boundary layer but 
there is evidence of jets much like those in the wake. 

Correlations measured in grid turbulence have been found 
to be highly anisotropic and consistent with the presence of vortex 
pair eddies. When a plane strain was applied to grid turbulence, 
the effect on the correlations suggested the presence of a stress 
instability similar to that postulated for the wake. 


1. INTRODUCTION 

‘Townsend’s (1956) measurements of the ‘transverse double velocity 
correlation’ in the turbulent wake of a circular cylinder have made it clear 
that a certain amount of order exists in the largest eddies in this flow. Using 
the downstream component of turbulent velocity, he found that the transverse 
correlation was positive for all values of the separation in the direction 
normal to the plane of the wake but was negative for large values of the 
separation in the direction parallel to the plane of the wake. ‘These 
observations were interpreted in terms of a set of large eddies, rotating 
about axes in the direction of the mean flow and with planes of circulation 
inclined so that the mean shear acted to decrease the area of circulation 
and increase the energy of the eddy. It was postulated that the average 
eddy was in energy equilibrium, gaining energy from the mean flow at the 
same rate as it was lost through the working of the Reynolds stresses of the 
more disorganized turbulence. 

This interpretation described the available data from the wake plausibly 
enough but the almost complete lack of comparable experimental observa- 
tions in other turbulent shear flows suggested the need for a more detailed 
examination of the correlation. This paper describes measurements of 
nine components of the double velocity correlation in a turbulent wake 
and a turbulent boundary layer. The results are used to form a new 
physical model of the largest eddies in the flow. 


2. EXPERIMENTAL PROCEDURE 


2.1. Coordinate system 


Cartesian coordinates are appropriate to all of the flows that are being 
considered. In the wake, the z-axis is parallel to the axis of the cylinder, 
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which lies in a direction normal to the free stream. The x-axis is in the 
plane of symmetry in the direction of the free stream and wx is increasing 
in the direction of the mean flow. The y-direction is normal to the plane 
of symmetry. 

In the boundary layer, the y-axis is normal to the wall and the x-axis 
lies in the wall in such a way that the direction of the mean velocity is parallel 
to the (x, y)-plane. The z-direction is also in the wall and normal to the 
v-axis. In the case of the square mesh grid, the x-axis is in the direction of 
the mean velocity, the z-axis parallel to the trailing set of bars and the y-axis 
parallel to the leading set of bars. 

‘The components of mean velocity in the x-, y-, s-directions respectively 
will be U, V, W, and the components of velocity fluctuation will be u, v, w. 
It is often convenient to use subscripts to distinguish between the directions, 
and the subscripts 1, 2, 3, will refer to the x-, y-, z-directions respectively. 


2.2. The double velocity correlation 

The quantity that kas come to be known as the double velocity correlation 
in turbulence theory, is the covariance in statisticians’ terminology. It will 
be written 


R’,,(x,x+r) = u, u, 


where u; is the 7-component of turbulent velocity at a point P with position 
vector x, and w’ is the j-component of turbulent velocity at a point P’ with 
position vector x+r. ‘The present discussion will be confined to cases 
where 7 = }. 

It is easier to measure the correlation coefficient 
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which does not require absolute determination of the sensitivity of the hot 
wires, but the denominator of this quantity is, in general, a function of r. 
It is therefore preferable to replace the denominator by the intensity at 
the fixed probe. For brevity in what follows, the quantities like 
R’(X, V5 14,2) 13)/u2, where u? and (x,y) are understood to refer to the 
fixed probe, will simply be written in the form 
Ri(hy V2 13)s 

and will be referred to as ‘correlations’ 

It is sometimes difficult to think of an eddy as a clearly defined structure. 
In general we are forced to rely on the intuitive idea that the turbulence 
consists of superimposed motions of various length scales and to define 
an eddy as a region over which the turbulent velocity is correlated. It may 
then be assumed that the covariance between the velocities at two points 
is unaffected by motions which have length scales much smaller than the 
distance between the points r. 

We may expect that the maximum length scale associated with the 
turbulence will be of the order of the largest length scale of the mean velocity 
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variation so it will be assumed that the covariance — 0 as r > ow. ‘This 
includes the assumption that any periodic motions which may be 
introduced into the turbulence do not extend very far in space without 
damping or change of phase. At the largest values of r for which the 
covariance is non-zero, the covariance will be dominated by the effects of 
the largest eddies present and at any smaller value of r, the covariance will 
receive contributions from all motions with length scales greater than r. 
In the particular case where there is a set of eddies, an order of magnitude 
larger than any of the other turbulent motions which contain appreciable 
energy, these eddies alone will determine the covariance over a considerable 
range of r at the end of the curve. 

The most obvious conclusions that can be drawn from a knowledge 
of the correlations concerns the backflow required by continuity. It may 
be shown (Townsend 1956, p. 14) that for two points in a plane, the 
correlation between the turbulent velocity components normal to the plane 
must be negative for some values of the separation vector r. For such 
values of r, the instantaneous velocity at the two points is, on the average, 
in opposite directions and if |r| is large, the correlation is dependent only 
on the largest eddies. We can thus gain some qualitative information 
about the large eddies quite easily, but even the whole correlation tensor 
is not by any means a complete specification of the motion. 

‘The one-dimensional spectra measured in the cylinder wake by Townsend 
(1950) show an appreciable peak at small wave-numbers fairly well separated 
on the wave-number axis from the broad peak which contains most of the 
energy inthe spectrum. It will be shown later that the observed correlations 
show considerable variations of character at large r. hese observations 
suggest the possibility of a fairly well-ordered set of large eddies with 
characteristic length scale large compared with the largest of the ‘energy 
containing’ eddies. 

Since we cannot predict a unique set of large eddies from the corre- 
lations, we may only hope to make a convincing argument for a rather simple 
set of eddies which will yield quantitative agreement with a number of 
observed correlations with an equal or smaller number of adjustable 
constants. Even if we can find such a fit, with a small number of constants, 
it is desirable to have further evidence for the chosen form. This might 
be a physical argument showing how the form could arise and should 
include at least qualitative agreement with all measured correlations. 

If the typical eddy can be expressed analytically, it is usually quite easy 
to calculate its contribution to the correlation. ‘The essential assumptions 
are as follows. 

1. ‘(he eddies are of a size much larger than any other component of 

the turbulence. 

2. ‘Vhe eddies are placed at random in the (x,)-plane. (In most 

shear flows, there will have to be restrictions on the location of the 
eddies in the y-direction. In a two-dimensional flow, the largest 


eddies may be comparable in size with the thickness of the turbulent 
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region so we must assume that all of the large eddies are located 
within a limited region on the y-axis.) 

3. The motions of adjacent or overlapping eddies are independent. 

(This ensures that any contribution to the correlation comes from 
a single eddy. ) 

Experience has shown that when a large set of measured correlations 
is available, only a very few kinds of simple eddy can yield correlations with 
the right sign at large r. ‘The sign of the expected correlation is easily 
determined without any calculation at all, simply by looking at a model 
of the assumed big eddies. One simply chooses two points on opposite 
sides of the model, separated in the direction corresponding to the 
correlation under consideration and notes whether or not the velocity 
component has the same direction at the two points. 


2.3. Experimental equipment 

The wind-tunnel used for measurement of wake and grid turbulence 
is a small closed return tunnel with a contraction ratio of 9:1 built and 
described by McPhail (1944). It has a working section 38 cm square and 
160 cm long. ‘The tunnel used for most of the boundary layer work had 
its origin in the duct used by ‘Townsend (1954) to distort grid turbulence 
by the application of a plane shear. It has undergone many changes since 
then but still retains the one-dimensional contraction and the motor section 
described by ‘Townsend. Between these, there now exists a duct 106 in. 
long and approximately 24 in. x 6 in. in cross-section. The cross-sectional 
area of the duct is constant for the first 50 in. and then increases 
gradually to compensate for the growth of the boundary layers. ‘The roof 
has been made adjustable in this latter region and has been arranged so 
that there is no longitudinal pressure gradient in the boundary layer on the 
floor. ‘There is a small trip fence on the floor at the exit of the contraction 
to ensure uniform transition and to thicken the boundary layer. 

The hot wires were made of platinum-cored Wollaston wire, etched 
after being soldered in place. ‘The diameter was 0-00025 cm (after etching) ; 
the length of the wire was typically 0-05 cm and never exceeded 0-1em. 
Two forms of sensitive element were used: a single straight wire placed 
normal to the stream (U-wire) for measuring the downstream component 
of turbulent velocity, and a pair of wires in the form of an X (X-wire) for 
the cross-stream components. Because of the need to bring X-wires into 
close proximity with each other and with a wall, it was necessary to construct 
several novel forms of wire holder. ‘The only special technique required, 
however, is patience. 


2.4. Accuracy of the correlation measurements 

The discussion that follows is confined to the part of the correlation 
curves for which ¢ is greater than 1-5 cm. As this is about 30 times the 
wire length, no wire length corrections have been necessary. ‘There is 
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also no need for precise compensation of each of the four hot wires. ‘The 
response of the wires is significantly non-linear in the boundary layer very 
near the wall but the effect on the correlation does not amount to more 
than about 5°.,, which is not important in the small correlations found in 
this situation at large values of r. 

The largest errors arise in the taking of averages. By averaging over 
a sufficient time, it is possible to determine the correlation at a given pair 
of points to a fraction of 1°,, but in order to save time in a long series of 
measurements, it is necessary to use averages which are sufficiently accurate 
for the purpose in hand but not necessarily the best that can be obtained. 
Generally the averaging time was adjusted to yield a probable error of about 
-Q-01 for correlations larger than 20°,, and decreasing for smaller corre- 
lations so that values of less than 0-01 have a probable error of about 0-002. 

Exceptions to the above estimates must be made for correlations with 
separations in the x-direction. In these cases, it is almost impossible to 
avoid some interference at the downstream wire from the wake of the 
upstream one. The procedure which has been adopted to minimize this 
effect has been to mount the sensitive part of the upstream wire on long 
lengths of Wollaston wire so that it stands out about 0-5 cm from the heavy 
supports. The track followed by the moving (downstream) wire passes 
within a few tenths of a millimetre on the side away from the supports. 
For values of r less than about 2 cm, the signal from the downstream wire 
shows some very high frequency components when viewed with an 
oscilloscope and peculiar effects are found in some of the measured 
correlations. ‘The correlations in this region are very sensitive to slight 
transverse movements of the downstream wire and are not very reliable. 
\t about r = 2 cm, the wake of the Wollaston wire appears to decay to the 
extent that its effects are negligible and the measured correlations become 
quite reproducible. ‘This condition lasts until about 6cm downstream 
where the wake of the wire supports in the leading probe begins to interfere. 
The errors in the observations are less serious 1n this region because this 
wake can be avoided by moving the trailing probe a short distance in a 
transverse direction. ‘The line joining the two probes changes direction 
by only two or three degrees. 


3. "THE WAKE OF A CIRCULAR CYLINDER 


3.1. Correlations in the cylinder wake 

The measurements to be considered have been made with the fixed 
probe at a distance of 533 diameters behind a circular cylinder at a Reynolds 
number of 1300 (U = 625 cm/sec, d=0-317cm). The mean velocity 
distribution was measured with a total head tube of 1 mm outside diameter 
and is given in figure 1. ‘The diameter of the cylinder (d) will be used as 
a unit of length and figure 1 serves to show the width of the wake in these 


units. The width of the intensity distributions can be compared with that 
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of the mean velocity distribution by reference to a comprehensive set of 
intensity measurements given by ‘Townsend (1956). 

The nine quantities R;(v;7,,72,73) are shown in figure 2 for cases 
where the separation r is in one of the coordinate directions. Each component 
has been determined with the fixed probe at three distances from the centre 
of the wake. One position is at the centre, one at 4 cylinder diameters 
from the centre which, at this value of x, is near the point of maximum shear, 
and one at 7-6 diameters which is near the outside edge of the wake but at 
a point where an appreciable signal can still be obtained. 
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Figure 1. The mean velocity profile in the cylinder wake at x/d S35 


These curves exhibit a wide variety of forms at large values of r. This 
suggests the presence of rather simple eddy structures because in a very 
confused turbulence the characteristics of any particular motion which 
does not appear frequently with a fixed orientation would be expected to 
disappear in the averaging, leading to smooth correlation curves with 
considerable symmetry. In isotropic turbulence, for example, continuity 
and symmetry require the ‘g(r)-type’ correlations (separation vector 
normal to the velocity component) to be identical and to be negative at 
some values of r. It is almost certain that these negative values would occur 
at the largest values of r for which the covariance is non-zero. The 
remaining correlations, which have the separation vector parallel to the 
velocity component, will be equal to each other. No proof has been given 
that they remain positive but this would seem to be the simplest physical 
situation and therefore perhaps the most likely. 

Examination of figure 2 will show that the observed correlations in the 
cylinder wake are very different from what is expected in isotropic 
turbulence. R;,(0,0,7) is negative for all values of y. When the 
fixed point is near the outside of the wake, R,,(0,7,0) has a minimum 
and R,,(0,r,0) goes strongly negative and then takes positive values as r 
increases. R,,(0,9,7r), Roo(r, 0,0), and R,,(r, 0,0) also have negative 
values followed by positive ones, at least for the larger values of y, but 
R,,(0,0,r) is independent of y and never negative. 
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3.2. The vortex pair eddies 

As a first attempt to interpret these curves, I searched for a set of a single 
kind of large eddy which would account for all of the observed correlations 
at large values of r. It is clear that the type of eddy chosen by ‘Townsend 
is not satisfactory because it is very long in the x-direction and its velocities 
are independent of x. This would lead to positive values of R,(r, 0,0) 
and R,,(r,0,0). There is a component of vorticity in the x-axis which 
would lead to negative values of R,.(0,0,7) which are not observed. It 
may also be noted that the correlations with separation in the y-direction 
all show significant non-zero values when the probes are in opposite sides 
of the wake. 

About the simplest kind of big eddy that we might hope to discover 
would be a body of fluid rotating about an axis, but not as a rigid body and 
with its circulation not necessarily in planes normal to the axis. ‘Townsend’s 
eddies are of this sort. The obvious difficulty in trying to produce the 
desired correlations from a distribution of these eddies is the observed 
values of R,.(0,0,7). There are certain configurations within the above 
definition of a simple eddy which can give negative values of R,;(0,0,7) 
but they require quite narrow restrictions on the orientation of the axis of 
the eddy and the planes of circulation with respect to the coordinate axes. 
It is doubtful if they could preduce negative values over a large enough 
range of r, and in any case these eddies seem inevitably to produce negative 
values of R,(0, 0,7). 

A more likely way to get negative values of R,(0,0,7) at large r, is to 
have pairs of simple eddies side by side and rotating in opposite directions. 
Again it is not possible to make a simple distribution of these eddies which 
will yield all the observed correlations at large r. Some progress can be 
made, however, if we assume that there exists a set of adjacent pairs of 
eddies in which the circulation is always in the (x,z)-plane. ‘These contain 
no v component of velocity and so only the affect R,, and R,, correlations. 
This type of eddy is illustrated schematically in figure 3. The arrows 
indicate the direction of the local velocity in the eddy and it can be seen 
that R,,(0, 0,7) will be negative and then positive as r increases; R,,(r, 0,0) 
will always be positive. ‘The bend in the axes of the eddy is needed for a 
qualitative explanation of the correlations with separations in the y-direction. 
It would be expected to produce a minimum in R,,(0,7,0) when the fixed 
point is in the outer part of the wake, and a negative region in R;,(0, 1,0). 

This eddy contains no v component and the correlations with separation 
in the (x,2)-plane will be independent of assumptions about the shape of 
the bend caused by the mean velocity if this is produced by a simple 
translation. ‘Thus we may calculate the correlations with separations in 
this plane from a model in which the axes of rotation are straight and in 
the y-direction. Such an eddy may be described by 
| 2 


u = A(1—2?2?)exp{ — 3a?(x? + 2?)} 


and w = Aa®xz exp{ — $a?(x*+27)}. 
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Assuming a uniform distribution of these eddies over the wake, we may 
calculate their contribution to the correlations. The results are as follows: 

R,,(r, 9,0) = a, exp( —4a?r?), 

Ry, (0, 0,7) = a, (1 — or? + Satyr exp(— far"), 

R33(r, 0,0) = R33(0, 0, r) = as {1 — 3a?r?}exp( — 4ar?). 
In these expressions, a, and a, are the proportion of u? and w? respectively 
which are contributed by the big eddies. If we set a, = a; = 0-13, and 
ad = ():25, these expressions fit the observations remarkably well as shown 
in figure 4. The order of magnitude of the constants is quite believable 


but it is a little hard to justify putting a, = a;. The intensities due to the 
big edilies are 





= — , — — 7 
= At 4 qw- = A2 
u;= A 2 and w, = ATS 
i 





Figure 3. Schematic model of the vortex pair eddy. 


We have then 

Now from values of v?/U? and w?/U? published by Townsend (1956, 
pp. 142-3), w/u2 is found to be 0-88 at y/d=4, but his X-wire was 
calibrated on the assumption that grid turbulence is isotropic and the 
most probable value is about 75°, of this (Grant & Nisbet 1957) which 


gives w?/u2 = 0:66. We now have 3 x 0:66(2 








:/wo*) = us/u®, or as = 3a). 
This is a rather uncertain calculation but it indicates that the assumed eddy 
does not fit quite as well as is implied by figure 4. However, even if the 
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Figure 5. Contour plot of R3,;(r,, 0, r3) expressed as percentage. y/d = 4. 
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condition a, = }a, is imposed upon the calculated correlation, there will 
still be quite a reasonable fit. 

It can be seen from figure 3 that although R,,(7,0,0) and R,,(0, 0,7) 
are negative for some values of r, Rj; should be positive for all values of 
separation in directions approximately at 45° to the x- and z-axes in the 
(x,z)-plane. ‘This has been tested by fixing one probe in the plane v/d = 4 
and moving the other over that plane. The result is shown in figure 5 
which is strong evidence that the eddy that has been assumed to be 
responsible for the correlations in the (x,z)-plane is qualitatively correct. 

The correlations with separation in the y-direction are more difficult 
to explain quantitatively by this model. ‘The search for an analytical 
description of the model is hampered by the need to make assumptions 
about the way in which the energy falls off with y, and about the shape of 
the axis of the eddy in the (x,yv)-plane. Since the wake is not even 
approximately homogeneous in the y-direction, we may also have to assume 
a frequency distribution describing the probability of the centre of the 
eddy being located at a given value of y. It cannot be assumed that all 
eddies are symmetrical about the central plane of the wake because that 
would require the correlation caused by the big eddies to be symmetrical 
about y = 0 also. For the fixed wire at y/d = 4, we would expect that 
R,,(0,7,0) and R,(0,7,0) would be approximately 0-13 at r/d= —8. 
This is clearly not so (figure 2). An attempt was made to fit the observations 
by assuming the eddies to be of the form 

u = A(1—a?2?)exp{— }a2[(x— by”)? +27] — Ja3(y—yo)*}, 

w= Aa?(x—by*)z exp{ — $a?[(x— by”)? +27] — 3a3(y —yo)*?. 
Here we have assumed a parabolic form for the bend in the eddy and an 
exponential decay of energy in the y-direction. In taking averages, yy was 
assumed to be uniformly distributed which is unrealistic in the outer part 
of the wake because, at y/d = 8, the total intensity is less than 13°, of the 
intensity at the centre while this assumption makes the intensity of the big 
eddies constant. ‘lhe result could not be made to give satisfactory agreement 
with the observations even for y/d = 4. Since we already have two constants, 


%, and 6, which are adjustable over a considerable range, there is not 


likely to be much profit in trying to refine the specification, as this would 
almost certainly introduce more unknown constants and_ eventually 
lead to agreement with the observations regardless of the real physical 
situation. 

R,, was measured with one probe fixed at the outside edge of the wake 
(v/d = 7-6) and the other moving over the (x, y)-plane. ‘The contours of 
zero correlation are shown in figure 6. ‘The shape of the region of positive 
correlation is about what is expected from the postulated model and 
indicates that the axes of the eddy are at about 45° to the coordinate axes. 
‘The turning points in the zero correlation contours are not at the centre 
of the wake but the reason for this can be seen if the data are plotted 
differently. 
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In figure 7 and several subsequent figures, the data are presented in 
a manner which is more useful for some purposes than the contour plot. 
The rmeasurements were made by fixing one probe at a point in the 
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Figure 6. Zero correlation contours for R3;(7,, 72, 0) with the fixed probe at y/d = 7-6 
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The scale of y/d refers to the moving probe. 
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Figure 7. Sections of the surface representing R33(7;, 72, 0) with the fixed probe at 
y/d = 7:6. The fixed probe and the straight lines followed by the moving 
probe are shown in their correct positions in the (7, 72)-plane. The correla- 
tion is plotted with these lines as axes. 


(x, y)-plane and moving the other along lines y = constant. In the figure, 
the position of the fixed probe, and the lines followed by the moving probe 
are shown in the (r,,7r2)-plane. Each traverse of the moving probe along 
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a line y = constant, determines a section of the correlation surface and 
these sections are plotted, using the corresponding lines as axes. ‘The 
scale for the correlation curves is shown at the left side of the figure. 

In figure 7, it can be seen that, at large values of r,, the curves have 
the shape which would be expected from a superposition of the correlation 
due to the vortex pair eddies and a negative correlation in the region around 
r, =. This is clearly the reason why R,,(0,7, 0) does not rise to something 
like 0-13 when the probes are at y/d = 8 and y/d = —8. The physical 
interpretation of this additional correlation is open to question, but it is most 
likely indicative of a tilt in the stress releasing motions which will be discussed 
in §3.5. The effect appears to be too localized to be caused by disorganized 
turbulence of this scale but the possibility of large eddies rotating about 
axes in the x-direction at the centre of the wake could only be ruled out by 
a detailed study of R,,(0,7.,r,) which has not been made. There is, 
however, no further positive evidence to suggest such a motion. 
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Figure 8. Sections of the surface representing R,,(7,, “2, 0) with the fixed probe at 
y/d = 7°6. 


The distribution ot R,, in the (x,y)-plane with the fixed probe again 
at y/d = 7-6 is shown in figure 8. The pattern is more complex than would 
be expected from the big eddies under discussion. It does, however, show 
a maximum correlation which is displaced in the upstream direction at 
the centre of the wake by about the same distance as the maximum 
correlation in the R,, plot. When the moving probe is in the other side 
of the wake, there is again a negative correlation, centred on the line x = 0, 
superimposed on the correlation due to the vortex pair eddies. Before we 
can attempt to interpret this, we must consider the correlations involving 
the v component of velocity. 


3.3. The mixing jets 
The R,, correlations for separations in the axial directions are shown 
in figure 2. The interesting points include the large negative values of 
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R39(r, 0, 0) | followed by positive values, the lack of negative values of 


R,,(0,0,r) and the fact that R,.(0,0,7r) extends for a comparatively short 
distance on the r-axis and is independent of y. Figures 9 and 10 illustrate 
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Figure 9. Contour plot of R,.(r,, 0, 73) at y/d = 4. 
%s 
10 Fined Moving Probe 
r / a\\ 
omit. / _——__ ],, 


4 
a 


0-8} 


| 
' 
y 
( 
\ 





— 





aA 
q 
ml Ps 
/ 
| 























a a — ‘ 
O2+ ss 
TT eee — 7-6 
i ET "ttt" 
O};- | 
~24 16 =e ° 7 6 
Va 


Figure 10. Sections of the surface of R»»(r;, 72, 0) with the fixed probe at y/d = 7-6. 
the distribution of R,, in the (x, v)-plane at y/d = 4 and in the (x, z)-plane 
with the fixed probe at y/d = 7-6. ‘The first of these figures clearly indicates 
rotation about axes in the z-direction. ‘The second suggests a structure 
L2 
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which is nearly periodic in the x-direction. In this case, the correlation 
is very asymmetrical when the moving probe is at the centre of the wake 
but there is not any shift in phase such as would be expected if a simple 
structure were being distorted by the mean velocity distribution. The 
periodic part, then, must arise locally or move along with the speed of the 
free stream. ‘The largest negative values of the correlation are centred 
roughly on the line along which R,, and R,, are most positive. This is 
probably only a coincidence as the positive peaks are also higher on the 
upstream side of the fixed probe and there is no obvious way to alter the 
model of the vortex pair eddies to include negative values of R,. along this 
line. On the whole, the R,, correlations appear to be the result of a different 
type of large scale motion with velocities mainly in the y-direction and with 
an approximate periodicity for short distances in the x-direction. 
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Figure 12. Sketch of a jet system in a late stage of development. 


At this point it is convenient to refer to experiments in which a circular 
cylinder emitting a small quantity of coloured fluid was moved through 
still water. The cylinder was } in. in diameter and 6in. long and was 
moved through a long glass-sided tank with cross-sectional dimensions 
6in.x6in. The ‘Reynolds number was about 800. The observer of this 
wake is in a position equivalent to that of moving with the stream in the 
wind-tunnel situation and he is able to see quite well the development of 
the outer part of the wake with time. Looking in the z-direction, the most 
obvious large scale motion in the outer part of the wake is a series of more 
or less regularly spaced ‘jets’ of turbulent fluid proceeding outward from 
the central plane of the wake. ‘Typical photographs are shown in figure 11 
(plate 1). From direct observation of the developing wake and from 
motion pictures, it appears that a series of jets in a late stage of development 
is part of a system similar to the one sketched in figure 12. The dashed lines 
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represent irrotational motions which rarely make a complete revolution. 
The structure is only approximately periodic, and as time goes on and 
successive ‘jets’ arise and decay, the typical ‘wavelength’ increases. 
This is obviously the sort of structure needed to produce the heavily 
damped periodicity in R,, but we must also consider the meaning of the 
asymmetry. It is observed in the tank that, when a jet is fully developed, 
the region of the turbulent core from which it arose has been displaced 
appreciably in the direction toward the cylinder. It may therefore be 
expected that the motion in the outer part of a jet is more closely connected 
to motions in the centre of the wake which are slightly upstream. The 
periodicity clearly indicates that the correlations extend over more than 
one jet, but if the structure of the jets is more regular than the distance 
between them, then any correlation which is produced mainly by the 
distribution of velocities in a single jet will be less ‘smeared out’ than one 
depending on having the probes in different jets. ‘The asymmetry of the 
negative regions is thus quite understandable and these general comments 
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Figure 13. Sections of the surface of Ro2(r,, 72, 0) with the fixed probe at the centre 
of the wake. 


would seem to apply to the outer positive regions as well but I am not able 
to offer any explanation of the asymmetry of the positive regions in terms 
of the detailed structure. Figure 13 shows the R,, correlation with the fixed 
wire at the centre of the wake. The positive peaks become asymmetrical 
when the movable wire is only a short distance from the wake centre, which 
indicates that they cannot be accounted for by the inward velocity at B 
(figure 12) which is a rolling-up of the end of the tongue of turbulent fluid 
and does not extend into the core of the wake. 

This model also suggests an explanation for the minimum in R,, along 
the line x = constant which passes through the fixed wire (figure 8). The 
large scale irrotational motions whose shapes are indicated roughly by 
dashed lines in figure 12 will contribute negative values of R,, with the 
fixed wire at the outside of the wake and the other at the wake centre. When 
the wires are at opposite sides of the wake and at the same value of x, then 
when one wire is at a place where the contribution to u from the large 
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irrotational motion is a maximum, the other is likely to be in a place where 
it is nearly zero. Upstream and downstream of this last wire, positive 
contributions to u will be found. ‘This situation would lead to negative 
values of R,, at this point but a positive contribution to the correlation is 
expected from the vortex pair eddies so that the observed values of R,, 
are qualitatively consistent with the assumed structure. 

The correlation R,. has a positive peak along the line x = constant, 
z = constant, which passes through the fixed wire (figure 10). For values 
of r, less than half the wake width, this may be attributed to small scale 
disorganized eddies but when the two wires are on opposite sides of the 
wake, it requires that the location of the jets on opposite sides of the wake 
should be correlated but sufficiently out of phase to produce positive values 
of Ry. The jets in figure 12 show the relationship that is envisaged and is 
commonly observed in the water tank experiments. 


3.4. The origin and maintenance of the vortex pair eddies 

The correlation R,;(0,0,7) has been selected as the single correlation 
that is most representative of the vortex pair eddies, i.e. the one that seems 
most unlikely to arise from some other kind of eddy. ‘This has been measured 
for values of y near the point of maximum shear at 250 and 533 diameters 
from the }-in. cylinder and at 546 and 1075 diameters from the +-in. 
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Figure 14. R,,(0, 0, r) at different distances from the cylinder. 
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cylinder. In all cases, U=625 cm/sec. The results are plotted in 
figure 14 and indicate a self-preserving increase of scale after 500 diameters 
at about the same rate as the increase of the mean velocity scale. 

It has been assumed that these eddies contain no v component of 
velocity. This is supported by evidence from the nine covariance components 
(figure 2). R,o(0,0,r) is always positive so that any v component other 
than a translational velocity of the whole eddy must come through a tilt 
of the axes of rotation in the (y,z)-plane. Even if this occurs, it is not of 
much interest since a tilt in that direction does not put the eddy in a position 
to gain energy through stretching by the mean shear. It appears likely 
that these eddies have been initially formed in the neighbourhood of the 
cylinder and are merely decaying with time. A determination of Ry, in 
the (x,z)-plane has been made with the fixed wire at a distance of 2 in. 
(16 diameters) behind the }-in. cylinder. The result is shown in figure 15. 
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Figure 15. Contours of zero correlation for R3,(r,, 0, 73) with the fixed probe 
at x/d 16. Correlations expressed as percentages. 


This plot has many of the characteristics of the corresponding one at 533 
diameters (figure 6). Here, however, there is a periodicity in the x-direction, 
patches of negative and positive correlation alternating with the Karman 
street frequency. This probably represents a folding of the Karman eddies 
as they break up. 
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There is some evidence that, when these eddies are destroyed, the first 
thing that occurs is the development of waves in the long vortices. Roshko 
(1954) found a periodicity in the R,,(0,0,7) correlation in the Karman 
street region behind a cylinder at Reynolds numbers over 500, and Hamma, 
Long & Hegarty (1957) have described a similar structure in the transition 
of aboundary layer. They towed through still water a flat plate with a small 
cylinder attached spanwise about 30 cm from the leading edge. They found 
that long vortices were regularly shed from the cylinder and that, at a high 
enough Reynolds number, these vortices quickly developed periodic waves 
in the (x,2z)-direction. ‘The plane of the waves then rotated so that the 
downstream maxima were carried away from the plate and the stretching 
of the vortex lines in the velocity gradient resulted in an amplification of the 
periodic disturbance leading to ‘horse-shoe vortices’. If a similar thing 
happens in the wake of a single cylinder, one might expect to find that the 
break-up of the street would leave the remains of horse-shoe vortices 
oriented in assorted directions—especially if the rotation of the wavy 
vortex carries some of it into the opposite side of the wake. 

There are other conceivable mechanisms but any folding or faulting 
can produce pairs of eddies rotating in opposite directions and it is probable 
that any eddy whose plane of circulation is not very nearly the (x, z)-plane 
will be quickly stretched out by the mean flow gradient and destroyed. 
Those whose orientation is right may at this time go through the energy 
gaining stage proposed by ‘Townsend but eventually they will be destroyed. 
A pair of eddies suddenly thrust into position with their axes normal to 
the plane of the wake and planes of circulation in the (x,z)-plane may at 
this stage have its middle torn out by the mean velocity gradient, but 
later, as the gradients are reduced, it may re-establish itself across the wake. 
It thus seems possible that the eddies with which we are concerned are 
simply selected bits of the Karman street. 

I tried to verify this idea by producing a wake which began without 
periodicity but I was unable to make such a wake. The bodies whose wakes 
were investigated were flat plates with holes or roughness to produce a 
turbulent boundary layer on the body. In all cases a large periodic 
component was found in the turbulence behind the trailing edge and the 
correlation R,,(0,0,7) was found to be negative at large r. 

Having failed to produce a wake without vortex pair eddies, I investigated 
the persistence of these eddies in a flow with no transverse mean velocity 
gradient. A grid of parallel cylinders of diameter 0-052 in. and spacing 
0-30 in. was used for this purpose. It is well known that the mean velocity 
variation disappears within about 30 mesh-lengths of such a grid which 


measured at distances of 477 and 1010 cylinder diameters downstream 
of this grid and the result is shown in figure 16. The correlation is negative 
at large values of r and is self-preserving when plotted as a function of (xd)'*. 
The solid line on the figure is the self-preserving curve for R,,(0,0,7) in 
the wake of a single cylinder and it coincides with the grid observations 


in this case is a distance of 10in. ‘The correlation R,,(0,0,7) was 
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at large values of r. It is to be expected that there will be differences at 
smaller values of r because there is now a larger characteristic length in 
the problem (the mesh-length) and much of the energy of the original 
mean flow variation must have gone into turbulence with length scales of 
this order. ‘This suggests that this turbulence, which has no mean velocity 
variation, contains vortex pair eddies which are the same size and growing 
at the same rate as those in the wake ofa single circular cylinder. If that is 
so, the vortex pair eddies are probably formed near the cylinder and the 
mean velocity gradient must play no part in their development. The stream- 
lines of the eddies must then be parallel to the (x, z)-plane. 
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Figure 16. Comparison of R,,(0, 0, 7) in the wake of a parallel cylinder grid with 
that in the wake of a single cylinder. 


Unfortunately, this argument does not satisfactorily account for the 
shape of the eddy in the (x, y)-plane that is suggested by the values of R33, 
(figure 6). With respect to a point near the outer edge of the wake, the 
centre is only retarded by a distance of 3-3 cm at 500 diameters downstream 
from the cylinder. ‘The mean velocity distribution is known over most of 
the wake (Townsend 1956, pp. 140-2) and it can easily be shown that an 
eddy with its axis straight and in the y-direction at x/d = 30 should have 
its centre retarded by about 10cm at x/d = 500. ‘This might be partly 
accounted for by the presence of more disorganized turbulence of comparable 
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scale in the (x, y)-plane. The effect could hardly be great enough to account 
for the whole discrepancy, however, and it is not yet safe to conclude that 
the vortex pair eddies all originate near the cylinder, although it is hard to 
escape that conclusion in the case of the parallel cylinder grid. This question 
will be considered again after the details of the mixing jets have been 
discussed. 


3.5. The mechanism of the mixing jets 

When a non-turbulent viscous fluid is subjected to a strain there are 
stresses set up in the fluid which are in many ways similar to those in an 
elastic solid but with the difference that the stress is proportional not to 
the strain but to the rate of strain. ‘There is, therefore, no energy of 
deformation and when the rate of strain becomes zero there is no tendency 
for the fluid to return to its original configuration. ‘The work done on the 
fluid has gone into heat and perhaps a change of pressure. 

When a turbulent fluid is strained, the situation is somewhat different. 
The viscous forces extract energy directly from the mean flow in the same 
way but there is an additional stress resulting from the presence of the 
turbulence. This has commonly been considered to have the same character 
as a viscous stress and it is accounted for in the mixing length theories by 
defining an additional viscosity known as the ‘eddy viscosity’. This 
procedure has been recognized to be physically unsound (e.g. ‘Townsend 
1954, Stewart 1956) although in spite of that it is sometimes found 
convenient. For the present purpose, the important difference between 
the viscous stress and the turbulent stress is that work done against the 
latter goes not into heat but into the production of turbulent energy and 
that an important part of the process is a redistrbution of energy among 
the turbulent velocity components which results in something very like the 
energy of deformation in an elastic solid. 

In general, the process by which the turbulence gains energy at the 
expense of the mean velocity gradient is not well understood. Batchelor 
& Proudman (1954) have made a theoretical analysis of the effect on a 
turbulent field of an instantaneous finite plane strain and they have given 
expressions for change of energy of the different components of turbulent 
velocity. ‘The first-order effect is simply a redistribution of energy among 
the turbulent velocity components to provide an increased intensity of the 
component in the direction of the compression and a decreased intensity 
of the component in the direction of the extension. ‘The second-order 
terms show an increase in the total energy which is shared, although not 
equally, by all the velocity components. This is an extremely simple 
situation and the effect of a finite and prolonged rate of strain on turbulence 
which is neither isotropic nor homogeneous is bound to be more complicated. 
It is to be expected, however, that there will be a large effect similar to the 
redistribution of energy described above as it is essentially an increase of 
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velocities associated with vorticity aligned in the direction of stretch and 
a decrease of velocities associated with vorticity in the direction of 
compression. 

In a plane wake, the nature of the strain is a combination of a rotation 
and a plane strain with principal axes of strain in directions at 45° to the 
usual coordinates. The resultant anisotropy of turbulent intensity is also 
greatest for axes aligned approximately in these directions. The subscript s 
will be used to indicate velocities and directions in a coordinate system based 
on the principal axes of strain, the relation to the usual coordinates is 
indicated in figure 17 for one-half of the wake. The relation between the 
velocities is given by 
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Figure 17. Relation between the ordinary coordinate system and the one aligned 
with the principle directions of strain for one half of the wake. 


The degree of anisotropy can conveniently be expressed by the quantity 
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which he calls the shear coefficient. He finds that in the fully turbulent 
part of the wake, not too near the centre, the shear coefficient is nearly 
independent of x and y and has a value of about 0-4. ‘To explain this, he 
postulated a ‘structural equilibrium’ in which an increase of anisotropy 
was opposed by a relative change in the rates of transport of energy along 
the spectra of the two velocity components. ‘This proposal has not been 
entirely satisfactory, particularly in the other situation where it has been 
used (see $5.2) and a different explanation of the constancy of the shear 
coefficient is attempted below. 
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When the equations of mean motion are expressed with respect to 
the principal axes of stress, only the normal components of the Reynolds 
stress are non-zero and for a shear flow in a steady state we have 
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Stewart (1956) has pointed out that the quantities péu2/ax, and 
pov* cy, have the same effect as pressure gradients and that pu and pv= are 
equivalent to the components of a pressure which in anisotropic turbulence 
are different in different directions. If the turbulence is homogeneous, 
these pressures are uniform and such a situation could be stable. On the 
other hand, an isolated blob of anisotropic turbulence in an irrotational 
fluid is clearly unstable. ‘The pressures on its boundary will be large where 
the boundary is normal to the direction of high intensity and small where the 
boundary is normal to the direction of low intensity. The irrotational 
fluid can only respond with accelerations which will result in an extension 
of the blob in the direction of high intensity and a contraction in the 
direction of low intensity, which is the sense required to equalize the 
intensities. 

The blob of anisotropic turbulence can be said to contain a stress acting 
in directions which tend to return it to the shape which it had before being 
strained if that was how the anisotropy was produced. ‘There is a clear 
resemblance to the stress in an elastic solid but there are two important 
ditferences. One is that the stress is only meaningful if it involves averages 
over volumes large compared with the scale of the turbulence. ‘The other 
is that to the first order, the energy available in the stress does not represent 
work done by the straining agent. It has arisen by a redistribution of 
turbulent energy and so may not bear much relation to the work dene in 
straining the fluid. 

‘The turbulence in a wake is under a continuous applied strain which is 
non-uniform and of different sign in opposite sides of the wake and it is 
bounded on both sides by irrotational fluid which is not being strained. 
‘That the anisotropy does not increase indefinitely is shown by ‘Townsend’s 
measurements of the shear coefficient, so the turbulent stress must be 
released by some means. ‘This could conceivably come about through 
an unequal rate of transfer of energy along the spectrum in the two velocity 
components or through the medium of secondary flows similar to those 
expected about the isolated blob of strained turbulence which was discussed 
above. It is probable that both mechanisms operate to some extent and the 
mixing jets seem to be the result of stress relief by secondary flows. 


If the bounding surfaces of the wake were plane on a scale an order of 
magnitude larger than the turbulence, the stress might be in equilibrium 
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because the resultant force would be everywhere in the same direction 
across the boundary. ‘The situation would probably be unstable, however, 
and we can make a guess about the physical nature of the flow resulting 
from the instability. If there is a bump on the bounding surface, and the 
turbulence in the wake has been uniformly strained, the pressure acting 
across the surface on one side of the bump will be different from that on 
the other. ‘The irrotational fluid can only see this as a high pressure region 
and a low pressure region and fluid will begin to flow from one side of 
the bump to the other. The type of flow that is envisaged is illustrated in 
figure 18(a). This might be expected to lead to a jet similar to those 
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Figure 18. The effects of distortion of a boundary between stressed and unstressed 


fluid. 


which have been observed. It might appear that a hollow in the turbulent 
fluid would induce a jet which bends in the opposite direction, but this is 
probably not so. Consider figure 18 (5). If a flow is set up from B to A 
in the irrotational fluid so that the flow in the turbulent fluid is simply a 
reflection of the turbulent jet, i.e. a jet of laminar fluid entering at A, then 
the fluid in the region a would have to expand in the direction a—c, which 
would always be up the gradient of normal turbulent pressure because 
an irrotational inflow from a and an outflow to 6 would serve to increase 
the stress at c. It would thus seem more likely that the expansion at a 
would be in the outward direction and contraction at 5 from the outer 
side. This would lead to a boundary like the dashed line in figure 18 (c) 
before the stress had been much reduced. We now have two bumps and 
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flow in the irrotational fluid from A to C or from A to B would be in a 
direction to relieve the stress deep in the fluid. 

Even if the flow in the outer fluid in figure 18(c) is from A to B, there 
is still an inclined surface to the right of B which could be the source of 
another jet and the jet arising from a single bump will also quickly cause 
a similar depression of the surface of the in-going side. It would seem that 
a roughly periodic set of jets might be expected to occur when a single one 
is started and this periodicity can be seen both in the photographs and the 
correlations (figures 10, 11). In the water tank it was common to see a 
set of from three to five jets on one side of a wake, all of approximately the 
same size and arranged approximately periodically. At a later stage, a 
new and larger set would arise in the same region of the wake. 

‘The boundary of a wake is never as simple as this model but it should 
be noted that in a fully developed wake, the boundary that matters is not 
the boundary of turbulent fluid but the boundary of stressed fluid. 
Although the outer part of a wake is very ragged after a series of jets has 
spent itself, the boundary of newly-stressed fluid may be quite regular, 
although diffuse, and some time may elapse before a suitable distortion 
of the ‘surface’ occurs to set off another series of jets. 

As the fluid is strained, the jets probably do not appear at any particular 
value of the stress, their initiation depending more on a chance configuration 
of the boundary. From observation of the wakes in water, I have formed 
the opinion that there is an irregular variation of the velocities and extent 
of travel of successive sets of jets. ‘The scale, however, appears to increase 
more regularly. Very soon after a jet begins, its effects probably reach 
right in to the centre of the wake and it is this length scale, the half-width 
of the mean velocity distribution, which probably determines the scale 
of the jets. The variation of Rj. over the (x, y)-plane was determined at 
x/d = 274. ‘The ratio of the ‘wavelength’ to that at x/d = 533 (figure 10) 
was found to be the same as the ratio of the widths of the mean velocity 
distribution at these two points (Townsend 1956, fig. 7.1). 

If we consider both sides of the wake together, it can be seen that their 
stress relieving action is not entirely independent. Referring to figure 19, 
a jet affecting the region A may draw fluid from the region B in the opposite 
side of the wake. ‘This would be expected to disturb the boundary of 
stressed fluid on the B side of the wake and produce a jet there as shown. 
It has been observed in the water tank that jets on one side of the wake are 
commonly accompanied by jets on the other side with their roots displaced 
slightly in the x-direction. 

The equilibrium value of the shear stress is probably an average over 
all stages of the build-up and release of stress. Why it takes the particular 
value it does, cannot be fully understood without an analytical treatment 
of the instability, but it probably depends on the average time required 
for a suitable perturbation of the boundary to occur. This model requires 
that the local value of the shear stress should vary in time and space, taking 
values ranging from near zero to something larger than the measured 
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average value. Evidence for this might be found in the probability 
distribution of the turbulent velocity components in the directions of the 
principal axes of stress. ‘Townsend (1956, pp. 153-4) has obtained values 
of the flatness factor of the distribution of u, and v,. Ata point away from 
the centre of the wake but not so far out that the intermittency is important, 
he finds that the flatness factors of u, and v, are 2-6 and 3-6 respectively. 
In the same region the flatness factors of u, v and w are very nearly 3-0 
which is the value for a normal distribution (‘Townsend 1948). 
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Figure 19. The relation between jets on opposite sides of the wake. 


This information by itself is not very useful but it does indicate that the 
distribution does not consist of patches of stressed and unstressed fluid. 
If that were the case, higher values of the flatness factor would be expected. 
The observed values are not inconsistent with the idea that the frequency 
distribution of velocities is derived from a field which contains a wide but 
continuous range of values of the stress. 


3.6. The rate of strain 

As the mean velocity distribution shows, there is a finite mean rate of 
strain over much of the wake and it is not changing very rapidly with distance 
from the cylinder. ‘The mechanism for releasing the excess stress must 
therefore operate without removing all the strain and without transferring 
it all to the outside of the wake. There are at least two reasons why the total 
strain is not removed with the accumulated stress. One of these is that the 
rotational part of the applied strain serves to rotate the principal directions 
of stress with respect to the principal directions of strain. This rotation 
might be as much as 15°. ‘The other reason is that any stressed fluid carried 
into the opposite side of the wake becomes fluid of negative stress in its 
new surroundings. ‘This fluid can then suffer some positive strain without 
building up any positive stress. We must aiso admit the possibility that 
some stress disappears through unequal rates of transfer of energy along the 
spectra of the different velocity components. It is not possible to estimate 
the magnitude of these effects quantitatively because the geometry of the 
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stress releasing motions is not known well enough. ‘There must certainly 
be some loss of total strain the in process so that, within a series of mixing 
jets, the downstream velocity near the centre of the wake is greater than the 
local mean velocity and the velocity outside the plane of maximum shear 
is less than the local mean velocity. As the jets die, these velocities on both 
sides of the plane of maximum shear will, on the average, approach the mean 
for the local values of (x,y) because of the drag of the fluid on either side 
of the region occupied by the jets. When averaged over the cycle of build-up 
and release of stress, the rate of strain within the region occupied by the 
jets is less than that which is indicated by the mean velocity gradient and the 
effect will appear in the mean velocity distribution as a gradual reduction 
of the gradients with increasing x. At the exteme edge of the wake, the 
immediate effect is an increase of the mean velocity gradient which, on 
decay, will produce an increase in the width of the mean velocity distribution. 


3.7. Relation between the vortex pair eddies and the mixing jets 

The form of the vortex pair eddies has been quite well established but 
the evidence regarding their origin and their relation to the mixing jets is 
somewhat conflicting. ‘hey seem to exist in the wake of the parallel cylinder 
grid, which has no mean velocity gradient. It is hard to see how such a 
highly ordered structure could arise within the turbulence by the action 
of the non-linear terms of the Navier-Stokes equations, so it would seem 
that they must, in this case at least, have been formed near the cylinder. 
The objection to this is that in the wake of a single cylinder, where the 
structure can be seen more clearly, the centre of the eddy does not appear 
to be retarded enough to have experienced the mean shear for the 
whole time required to pass from near the cylinder to the point of 
observation. 

It is possible that the phenomenon occurring in the cylinder wake is 
different from that in the wake of the grid. On this supposition, one could 
argue that the vortex pair eddies arise in the wake in its fully developed state. 
The likely mechanism would be that the eddies form about a tongue of fast 
irrotational fluid which is drawn into the wake as a part of the mixing jet 
action. ‘There is, however, no sign of vortex pair eddies in the outer part 
of the turbulent boundary layer. 

There is not enough evidence to resolve the problem with certainty 
but I think it most likely that these eddies do arise in the region of the 
vortex street. [he reason for the small retardation of the centres of the 
eddies may be that they have a tendency to set off a mixing jet. If they 
have such a tendency, they will, on the average, be found superimposed on 
weak jets which appear before the stress vector has rotated much with 
respect to the strain. ‘They would then have less than the average mean 
velocity gradient and would also be more likely to survive in a coherent 
form than if they were subject to some of the very large jets which draw 


irrotational fluid almost into the centre of the water wakes. 











“si 


The large eddies of turbulent motion 17 


4+, "THE TURBULENT BOUNDARY LAYER 

4.1. Correlations in the boundary layer 

Although the study of the large scale structure of the wake has been 
described first, the measurements were actually made after a similar but 
less extensive series in a turbulent boundary layer. The boundary layer 
experiments were not designed in the light of the experience with the 
wake, which is a simpler flow, but can be profitably discussed that way. 

The wind-tunnel has been described briefly in $2.3. The observations 
were made at a distance of 80 in. from the trip fence, and the mean velocity 
profile at that point is given in figure 20. ‘The total thickness is about 
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Figure 20. Mean velocity profile in the turbulent boundary layer. 


5-5 cm and the thickness 5), which is the distance from the wall to the 
point at which U,—U=U_, is 3-83cm (U = mean velocity, U, = free 
stream velocity, and U2 = wall stress/fluid density). The wall stress has 
been determined from the slope of the logarithmic part of the curve on the 
assumption that it is represented by 


U = ss | tog ee +4] 


and that K = 0-41 (Townsend 1956, p. 242). This profile agrees very 
well with most published ‘flat plate’ boundary layer profiles when plotted 
in the form (U,—U)/U_ vs y/5). For example, it is in agreement over the 
whole layer with those collected by Townsend (1956, pp. 244-5). 

The nine correlations are given in figure 21. These curves differ 
considerably from the corresponding ones for the wake. Perhaps the 
most noticeable thing is the large decrease in most of the length scales as 
the fixed probe moves nearer the wall. Because of this it 1s best to consider 
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the correlations in the outer part of the layer separately from those near 
the wall. 


4.2. The observations in the outer part of the layer 

The observations in the outer part of the layer are collected together 
in figure 22. Most of them are taken with the fixed wire at y/d, = 0-66 
which is near enough to the wall for there to be not much intermittency. 
The curves approach zero over a considerable range of values of r, indicating 
a wide range of length seales for different combinations of velocity component 
and direction of separation in the largest eddies. The correlations for the 
w component are highly anisotropic even at small values of r. 
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Figure 22. The correlations in the outer part of the layer. 


The backflow corresponding to the different velocity components at 
the fixed probe can be located roughly. The curves for R,,(0,7r,0) and 
R,,(9,0,r) show the R,, correlation along two lines in the plane normal 
to the u component of velocity. It can be seen that the back-flow contributed 
by the largest eddies takes place mostly at values of r in the z-direction. 

M2 
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Similarly, for the v7 component it occurs mostly in the x-direction and for 
the w component in both the x- and y-directions although at different 
values of r. 

It is tempting to try to use these back-flows individually or in pairs 
to infer the nature of the large motions. Not much can be done, however, 
beyond statements that certain simple large circulations are not dominant. 
As in the wake, the trouble is that even the nine correlations are only a small 
part of the specifications of the motion and cannot be expected to yield a 
unique solution without the aid of further evidence or hypotheses. 

Although these correlations are quite different from the corresponding 
ones in the wake, their complication again indicates a considerable amount 
of order among the big eddies, and it is worth searching first for a set of 
large eddies of a single kind which will yield all the nine correlations at 
large ry. ‘Townsend’s large eddy structure will produce the observed 
values of R,,(0,7,0) and R,,(0,0,7) quite reasonably but since the 
directions of the velocities in these eddies are independent of x, they 
cannot yield negative values of R,.(r, 0,0) and R,,(r, 0,0) even if they are 
reduced to finite length. ‘They cannot therefore be used as a detailed 
picture of the large scale motion. 

There is no suggestion of vortex pair eddies, notably because R,,(0, 0, r) 
is clearly positive at all values of r. A number of other kinds of simple 
eddy were considered but there is not enough evidence to support any of 
them as a dominant structure. 

The possibility of mixing jets similar to those observed in the wake 
would seem a likely one and most of the signs are there. R,.(0,0,7) goes 
very slightly negative this time but the negative values are so small that 
the character of the curve bears a clear resemblance to the corresponding 
one in the wake. ,,(r, 0,0) takes negative values as large as 5”, but there 
is no sign of a periodicity here. R,,(0,7,0) is positive except for one point, 
out where the intensity is very low. R,.(0,7,0) is always distinctively 
positive as in the wake. 

All the relevant components are thus consistent with the mixing jet 
idea except for the lack of periodicity in R,,(7, 0,0). A measurement of R35, 
over the (x, y)-plane is given in figure 23. ‘The fixed wire was at a distance 5, 
from the wall and the other wire was moved along three lines between the 
tixed wire and the wall. ‘The sections of the correlation surface have an 
asymmetry very much like that of the corresponding correlation in the wake 
(figure 10). ‘here is no periodicity, but the curves have long flat negative 
tails which suggest motions with a wide range of sizes. 

The extremely large rates of shear and the presence of the wall produce 
turbulence of a very small scale in the inner part of the boundary layer and 
since this scale is very much smaller than the thickness of the layer, it is 
likely that the turbulence undergoes stress releasing motions in this region 
which, although large compared with the local turbulence scale, are still 
small compared with 6). Some of these will appear in the outer layer and 
the stress releasing motions will have quite a broad spectrum there. 
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Klebanoff (1955) has measured energy spectra in a boundary layer, and, 
except very far out in the intermittent region (v/d) = 1-5), he found no 
peak in the spectrum. ‘This is in contrast with the clearly defined two- 
component structure found by Townsend (1950) in the wake spectra. 
The asymmetry of figure 23 indicates that the individual motions have 
a character similar to those in the wake, i.e. they go in the same direction. 
The type of instability involved is about the same as in the wake but the 
part of the unstressed fluid in that model is often played in the boundary 
layer by fluid with smaller stress and this leads to a much more confused 
motion. 
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Figure 23. Sections of the surface of R22(7;, r2, 0) with the fixed probe at y = 8p. 


4.3. The motion near the wall 


Figure 24 shows the nine correlations at points near the wall but in the 
logarithmic region of the mean velocity curve. The distance from the wall 
varies between 0-0346, and 0-0765). A representative value would be 0-055. 
Comparison with figure 22 shows large differences both in scale and character. 

There is no support here for the suggestion that the dominant form of 
the energy containing motions is the ‘attached eddy’. ‘The very great 
extent of R,,(r,0,0) compared with R,,.(r,0,0) and R,3(r,0,0) does not 
seem consistent with a simple eddy as we have used the term. The very 
small extent of R,.(0, 0,7) strengthens this view. 

These observations are consistent with the presence of stress releasing 
motions similar to those proposed for the outer part of the layer. The 
scale of these motions would be small compared with those originating in 
the outer part of the layer but large compared with the scale of the 
turbulence that contains most of the energy locally. If a series of jets of 
this size occur, lined up in the direction of the stream, then, because of the 
large velocity gradient, there will be a large contribution to the u component 
of turbulent velocity. This contribution will be a function of y but will 
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be nearly independent of x over the whole length of the series of jets. 
Such a structure could produce positive values of R,,(r,0,0) at values 
of r,; many times the distance from the wall. 

This picture is supported by several of the other correlations. The 
back-flow corresponding to the uw component of velocity at the fixed point 
occurs only in the z-direction, at least as far asthe bigger scales are concerned. 
This implies stress relieving motions originating very near the wall, perhaps 
involving the boundary of the laminar sublayer. ‘The motion would be in 
the nature of an outward eruption originating near the wall. 
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Figure 24. The correlations with the fixed probe near the wall. 


The extremely short extent of R,.(0,0,r) means that the series of 
motions giving rise to the extended correlations in the x-direction is quite 
narrow in the z-direction, and its movement in the y-direction is not 
compensated by an opposite movement or even an extended drift in the 
region to the side of the occurrence. Rather, the back-flow appears to be 
in the x-direction, which would be expected in a stress releasing motion or 
a sequence of them lined up in the x-direction. 
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The negative values of R,,(0,0,7) would then represent flow toward 
the displaced fluid on both sides because, very near the wall, rotational 
motions involving the w component of velocity might take the place of the 
irrotational flow from the other side which was postulated in the wake. 
The back-flow for the w component occurs only in the y-direction, which 
again suggests that the displaced fluid is long in the x-direction. When 
the fixed probe is moved a short distance away from the wall, R;3(r, 0,0) 
does become negative (figure 21), indicating a loss of order in the large 
scales of the w component. Slightly further out, at about y = 0-256), 
R,,(0,0,7) ceases to be negative. Beyond this level, it may be assumed 
that the nature of what may be locally described as large scale motions, is 
essentially the same as that described for the wake but with stress releasing 
motions of assorted size originating over a considerable range of values of y. 

This description of the motion near the wall is supported by photographs 
published by Einstein & Li (1956), showing coloured fluid from the laminar 
sublayer erupting into the turbulent fluid over a limited region of the wall. 
Unfortunately, the present experiments do not yield any information 
relevant to their proposal of an intermittent instability in the sublayer. 
‘The motions proposed here could be caused by the rapid production of 
shear stress which would occur just after a sudden transition in the sublayer. 
On the other hand, they might originate in the fully turbulent fluid but 
be vigorous enough to stir up part of the sublayer. 


4.4. Space time correlations 

Some suggestions can now be made about the interpretation of the 
space-time correlations measured in a boundary layer by Favre, Gaviglio 
& Dumas (1957). ‘The first of the two most outstanding points among 
their results is that for two probes separated in the y-direction, the time 
delay for maximum correlation increases, as the probe pair is moved 
downstream, at a faster rate than can be accounted for by the mean velocity 
gradient. ‘The other is that, for one probe fixed and the other taking 
positions over the (x, y)-plane downstream, if values of r, and the time 
delay were adjusted for maximum correlation at each value of r,, the locus 
of the points of maximum correlation was a curve which moved away from 
the wall much faster than a mean streamline. 

The first of these observations is to be expected if the large scale motions 
in the boundary layer are essentially the same as those which have been 
found in the wake, except that there is nothing equivalent to the interaction 
between the two sides of the wake. ‘That is, they consist of ‘eddies’ always 
rotating in the same direction, which is such that on the outer side they are 
moving faster than the local mean velocity and on the wall side slower than 
the local mean velocity. ‘The second observation indicates that (as has 
been assumed here) the outgoing part of the motion is more coherent than 
the ingoing part so that it can, when it is convenient, be thought of as a jet. 
The line of maximum correlation in time and space is thus something like 
an average streamline for particles within jets. The velocity of such particles 
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will, of course, be expected to maintain a non-zero correlation for much 
longer periods than particles not in jets. 


5. GRID TURBULENCE 
5.1. The large eddies of grid turbulence 

The presence of vortex pair eddies in the cylinder wake suggested 
that something similar might be found in grid turbulence which, at its 
origin, consists of a combination of cylinder wakes. Much of the experimental 
work on grid turbulence has been done with the assumption that, at some 
distance from the grid, the motion is an accurate approximation to theoretical 
isotropic turbulence. In the theoretical case, all nine correlations are 
determined by a single scalar function f(r), which in the present notation is 
R,,(7, 0,0) = R,o(0,7r, 0) = R33(0,0,7). The other six non-zero correlations 
have a single form g(r), which is a function of f(r) (Batchelor 1953). When 
this is assumed to be true, complete experimental information on the 
correlation tensor is obtained by determination of f(r) in the form R,,(r, 0, 0), 
but it is a simple matter to measure g(r) directly using the downstream 
component of velocity and this has also been done (Stewart & ‘Townsend 
1951). It has been shown by Grant & Nisbet (1957) that the turbulence 
is not isotropic in the sense that 72 < u?. The idea of isotropy with regard 
to the correlation components was thus in question and a direct test was 
desirable. 

The coordinate directions which have been chosen for the grid are those 
of the trailing set of bars. If the vortex pair eddies are to be found here, 
they will presumably come from the wakes of each set, so a single correlation 
will be a combination of contributions from each set of vortex pair eddies. 
Table 1 gives the nine correlations in the grid notation, and the components 
which they include in the notation appropriate to the two different sets of 
wakes. Also given are the signs of these correlations at large r and the sign 
that might be expected in the grid turbulence. ‘The results of some 
measurements of R,,(0,7,0), Rs3(0,0,7), Ro(0,0,r) and R,3(0,7,0) are 
given in figure 25 for a l-in. grid at 80 mesh-lengths and a 2-in. grid at 
30 mesh-lengths. In every case, the sign of the correlation at large r is 
different from that which is expected in isotropic turbulence. It appears 
that the large scale motions retain some of the anisotropic characteristics 
of the original wakes for a very long time; probably throughout the life 
of the turbulence. 


5.2. The effect of a plane strain on the large motions of grid turbulence 


Some years ago, ‘l’ownsend (1954) made a study of the effect of a plane 
strain on grid turbulence. He constructed a duct of constant cross-sectional 
area which, when placed downstream of a grid, compressed the turbulent 
fluid in one transverse direction and extended it in the other. He found 
that, in terms of coordinates aligned with the principal axes of strain, the 
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structure parameter A, = (v?—w*)/(v?+w?) became constant after a total 
strain ratio of about 2, The value of K, at this point was about 0:4, which 
is very close to the equilibrium value of the shear coefficient in the turbulent 
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Figure 25. Correlations measured in grid turbulence. 


‘Townsend suggested that the cause of the equilibrium is that the crowding 
together of the favoured eddies causes an increased rate of transfer of energy 
from these to smaller eddies. As he pointed out, the difficulty with this 
explanation is that when the distortion ceases, the anisotropy decreases 
only very slowly. ‘This hypothesis would require that eddies of a given 
energy can be located relative to each other in such a way that the energy 
transfer is normal but if any energy is transferred from one component of 
turbulent velocity to the other, it is immediately passed on to more disordered 
turbulence. 

I reconstructed the wind-tunnel with the distorted duct but with a 
longer parallel flow section after the distortion (figure 26). The coordinates 
used are those of the trailing set of bars in the grid which is oriented so 
that y is in the direction of compression in the distorting section and z is 
in the direction of extension. ‘The origin of x is at the grid. In this tunnel, 
I was able to measure A, at a distance of 26 in. after the end of the distortion. 
The turbulence was not very homogeneous but an average of seven readings 
at different points over the tunnel section was K, = 0-397 with a standard 











The large eddies of turbulent motion 187 


deviation of 0-024. There is no doubt that the anisotropy is decreasing 
only very slowly (figure 27). 
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Figure 26. The wind tunnel with the distorting duct which applies a uniform plane 
strain to grid turbulence. 
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The similarity between the behaviour of K, in this flow and in the wake, 
suggests that the mechanism which maintains a constant average stress is 
similar in the two cases. It is hard to predict the form of the stress relief 
in detail, but it must occur in patches or cells and the typical motion would 
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be something like that described for a ‘blob’ of anisotropic turbulence in 
an irrotational fluid. If the turbulence were perfectly homogeneous, it 
might be in equilibrium even for large stresses, but in practice there is 
probably always enough variation of intensity to make an instability possible 
for some level of stress. 

One thing that can be expected is that the stress relieving motions will 
be in the (y,z)-plane and be considerably larger than the typical scale of 
the turbulence. ‘This should lead to negative values of R,.(0,0,r) and 
R3,(0,7r,0) at large values of r. The result of a measurement of R,,(0, 0,7) 
at a point three inches before the end of the distortion is given in figure 28. 


2 





Figure 28. R..(0, 0, 7) near the downstream end of the distorted duct. Fixed 
probe at x = 57in. 


This point was 154 mesh-lengths from the grid and, in the absence of a 
strain, we would expect a curve not much different from those shown in 
figure 25(b). ‘The scale has increased by more than a factor of two and 
there are no negative values. ‘The increase of scale in this component is 
caused by the extension of eddies with vorticity in the z-direction and the 
stress releasing motions probably have a similar scale so that their negative 
contribution is masked by the effect of the long eddies. 

The behaviour of R,,(0,7,0) is more interesting. Figure 29 shows 
this correlation at several points in and after the distortion. Before the 
distortion this component is also positive at all values of r but at x = 39 in., 
which is 19 in. after the beginning of the distortion, it is negative foie values 
of r greater than one mesh-length. The equilibrium structure has not yet 
been reached at this point. Eighteen inches further downstream, at 
v = 57in., the scale has decreased considerably, and at large r the correlation 
has become still more negative. That point is in the equilibrium region, 
three inches before the end of the distortion. At x = 75-5, which is 15-5in. 
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after the end of the distortion, the scale has increased again but there has 
not been much change in the correlations at large r._ Fourteen inches 
further downstream at x = 89-5 in., however, the negative values of the 
correlation have decreased distinctly. 

In very disorganized turbulence, e.g. isotropic turbulence, this 
correlation is expected to be negative and one of the reasons why it 
changes sign in the distortion may be the suppression of the vortex pair 
eddies that make it positive in ordinary grid turbulence. ‘This is probably 
not enough though, because the negative values are quite large and extend 
over a large range of r. The decrease in the negative values after the 
distortion ends is to be expected if the stress relieving motions initiated 
in the equilibrium region of the distortion are completing their lives 
within twenty or thirty inches after the end of the distortion and not many 
new ones are being formed. 
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Figure 29. R;,(0, 7, 0) in and after the distortion. 


The observations are thus consistent with the main predictions of the 
hypothesis of stress relieving motions but it is necessary to assume that 
when the stress is not too large it is stable, otherwise it is hard to account 
for the very slow decrease of the anisotropy after the distortion. In the 
discussion of the wake, it was suggested that the time required for an unstable 
situation to develop a mixing jet governed the observed average value of 
the stress. Here we do not have the large intensity gradient of the wake 
and it may be that when the stress is below a certain value it acts to make the 
turbulence more homogeneous and is thus stable. 

The flatness factors of the distributions of v and w vary considerably 
over the section of the tunnel but mean values at the end of the distortion 
are for v, 2:97, and for w, 3-09. ‘Twenty-six inches after the end of the 
distortion, the corresponding values are 2-93 and 3-26. The standard 
deviations were about 0-02. ‘These differ from 3 in the same directions 
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as the corresponding quantities in the wake. Again, I do not understand 
why one is larger than 3 and the other smaller, but they show that the 
turbulence does not consist of patches of stressed and unstressed fluid 
even after the distortion has ceased, which suggests that if the stress is not 
too large, it is stable to large stress relieving motions. 

It thus seems that the values of the shear coefficient in the wake and in 
distorted grid turbulence are established by different mechanisms and the 
fact that the coefficient takes the same value in the two cases is only a 
coincidence. A constant, but different, value of the shear coefficient is 
found in several other kinds of turbulent flow. In fully developed 
two-dimensional channel flow (Laufer 1950) and in the turbulent boundary 
layer (Klebanoff 1955), the value is nearer 0-3 than 0-4. For an axisymmetric 
contraction, Uberoi (1956) gives data from which the coefficient K, can 
be calculated and it is found to reach an apparent equilibrium value of about 
0-8 after a total strain ratio of about 2 in the (x, y)-plane. 


I am indebted to Dr A. A. ‘Townsend for helpful discussions during 
the course of this work, and to the Defence Research Board of Canada for 
financial assistance. 
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An approximate measurement of the ionization time 
behind shock waves in air 


By BRYAN NIBLETT and VERNON H. BLACKMAN* 


Palmer Physical Laboratory, Princeton University 
(Received 7 August 1957 and in revised form 9 January 1958) 


SUMMARY 

A hydromagnetic shock tube has been used to obtain an 
approximate measurement of the time to reach equilibrium 
ionization behind shock waves of Mach number 11 to 17, moving 
into air at a pressure of about 1 mm of mercury. This ionization 
time decreases with increasing Mach number. Experimental 
results are presented as a graph of the ionization time vs Mach 
number. ‘The principal source of error in the measurements is 
the attenuation of the shock and the results indicate a lower limit 
for the ionization time in air. 


INTRODUCTION 

In recent years there has been considerable interest in the rate at which 
equilibrium ionization is reached behind shock waves in gases. Petschek 
& Byron (1957) measured the ionization rate for argon at shock velocities 
from Mach number 10 to 18; their results indicate that the dominant 
ionization process is collisions between electrons and argon atoms, and 
that there is good agreement between experimental measurements and 
theoretical calculations of the ionization rate. ‘Turner (1956) has made 
similar measurements in xenon gas. In a theoretical paper, Bond (1957) 
also concludes that the electron—atom process is the dominant one. 

No measurement of the ionization time behind shocks in air have been 
reported so far despite the importance of this information in aerodynamic 
problems. Lamb & Lin (1957) measured the electrical conductivity of 
shock heated air, and concluded that, at the densities used, the time to 
reach equilibrium ionization is short compared with the time interval 
which can be resolved in their experiments. ‘The present note reports a 
rough measurement of the ionization time in air behind shocks in the 
Mach number range 11 to 17 and at an initial pressure of about 1 mm of 
mercury. 

EXPERIMENTAL METHOD 

Plane shocks in air were generated using a hydromagnetic shock tube. 
This apparatus has already been described by Blackman, Niblett & Schrank 
(1957). It consists of a one-inch diameter glass tube with a single-turn 
copper coil wrapped around it. A low-inductance capacitor (rated at 


* Now at Giannini Research Laboratory, Santa Ana, California. 
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1-6 microfarads and 0-025 microhenrys) is discharged through the coil 
and the rapidly varying magnetic field (70 kilogauss per microsecond) 
induces an electric field which breaks down the air and produces circular 
current loops. An analysis of the initial process has been described by 
Blackman & Niblett (1958). Shock waves, moving outward in both 
directions from the coil, are produced by the rapid expansion of the gas 
due to Joule heating and the interaction of the magnetic field with the 
current loops. ‘The motion of the shock wave in the horizontal tube is 
photographed using a drum camera; photographs taken with a vertical 
slit show that the shock becomes plane within about three diameters of the 
coil, and photographs taken with a horizontal slit and drum camera yield 
the shock velocity. A typical streak photograph is shown as figure 1 
(plate 1). With the capacitor initially charged to 20 kilovolts and the air 
at 1 mm pressure, the luminous front achieves its maximum velocity of 
about Mach 90 at 1 cm from the coil and decays, rapidly at first and then 
more slowly, to Mach 10 at about 30cm from the coil. Over the range 
of Mach numbers 11 to 17 the rate of decay of the front is about 0-5 Mach 
numbers per cm. 

The shock front itself is not visible on the photographs, but is followed 
at some distance by a region where the luminosity increases rapidly. It is 
assumed that the radiation is produced mainly by processes which involve 
free electrons. ‘This assumption is reasonable since in the range of 
photographic sensitivity the radiation is produced principally by oxygen 
free-free and oxygen free—bound transitions together with nitrogen band 
spectra (see Keck et al. 1957) which are likely to be excited by free electrons. 
The luminous front is thus taken to be the point at which ionization 
equilibrium is attained. For strong shocks in which the ionization time 
is very short the luminous front and the shock front are coincident but 
as the wave decays the shock front draws progressively further ahead of 
the luminous front. ‘The intermediate dark space is the region in which 
the gas ‘relaxes’ to its equilibrium condition. 

The vertical distance on the film (i.e. parallel to the time axis) between 
the shock front and the luminous front is taken as the measure of the 
ionization time. ‘The essential problem, therefore, is to determine on each 
film the locus of the shock front. The position of the shock front was 
determined by reflecting the wave from the plane face of a one-inch diameter 
metal slug placed in the tube at appropriate distances from the coil. When 
the shock wave is reflected from the face of the slug the enthalpy of the gas 
is approximately doubled and the ionization time decreased. Over the 
range of these experiments the ionization time is reduced by a factor of 
about five on reflection and hence an error of about 20°, is incurred by 
assuming the reflected shock and luminous front to be coincident. 

A series of shocks were fired into air at pressures between 1 and 2 mm 
of mercury with the reflecting slug placed at varying distances from the 
coil. Fresh air was used for each experiment since the discharge may 
change the composition of the gas in the tube. A typical photograph 
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Figure 1. Typical streak photograph of a shock produced in the hydromagnetic shock tube. Distance 1 


in the horizontal direction and time in the vertical direction. ‘The coil is on the left and the vertica 
lines are markers 3 cm apart. Maximum Mach number is about 90 and occurs 1 cm from the coil. 
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ionization time behind shock waves in air, Plate 2. 
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Figure 2. This photograph shows a shock reflected from a metal slug. At the reflecting surface the 


incident Mach number is 11:1, and p,7 = 0:56, in units of cm of mercury times microseconds. 
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is shown as figure 2 (plate 2); the reflected shock, appearing earlier in time 
than the incident luminous front, can be seen clearly. Ifa series of pictures 
is taken with the same gas pressure and the same capacitor voltage, the 
streak records can be accurately superimposed and a ‘master’ diagram 
constructed showing the locus of the luminous front and the shock front. 
Such a diagram is shown in figure 3. The Mach number of the shock can 
be obtained from the angle (as indicated in the figure) and the film speed; 
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Figure 3. The position of the shock and luminous fronts as a function of time. The 
diagram is used to measure 7 and the Mach numbers. The path that a gas 
particle follows is indicated by the dotted line. ‘The scatter of the crosses, 
which represent the position of the shock front, gives a measure of the experi- 
mental precision. 


and the distance along the time axis between the shock and luminous front 
is a measure of the ionization time. ‘This measured time +7 must be 
multiplied by the mean density ratio across the shock, which is about 
nine in these experiments, to give the true ionization time; i.e. the time 
for a given sample of gas to reach equilibrium ionization along a particle 
path. 
RESULTS 

The results are plotted in figure 4 in terms of log(p, 7) vs Mach number, 
where p, is pressure in cm of mercury and 7 is in microseconds. The 
points fall roughly on a straight line when plotted in this way. Each point 
represents a value of 7 and Mach number measured at the same distance 
from the coil. 

The most serious source of error is the attenuation of the incident shock. 
The gas sample which reaches equilibrium ionization at a distance x, from 
the coil (see figure 3) has passed through the shock further upstream at a 
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distance x, where the Mach number was greater. ‘The Mach numbers 
plotted in figure 4 are therefore too low. The effect of this error has been 
avoided as much as possible in the following way. Measurements of the 
ionization time in argon have been made using the present technique and 
compared with the more accurate results of Petschek & Byron. Agreement 
was found to be reasonably good for short ionization time but increasingly 
poor for long times. Accordingly, results for air are included in figure 4 
only for the range of ionization times over which there is reasonably good 
agreement between the two sets of argon measurements. 
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Figure 4. The experimental values of the ionization time in air plotted against Mach 
number. The ionization time is multiplied by p, since the ionization rate 1s 
closely proportional to the atom density. 


‘The present work was supported by the Office of Naval Research, 
Washington, D.C. One of us (B. N.) wishes to acknowledge his indebted- 
ness to the United Kingdom Atomic Energy Authority for a grant of leave 
making possible his collaboration in this investigation. 


REFERENCES 

BLACKMAN, V. H., Nisiett, G. B. F. & ScHRANK, G. 1957 Bull. Amer. Phys. Soc. 2, 
216 (abstract). 

BLACKMAN, V. H. & Nisiett, G. B. F. 1958 Proc. 2nd. Lockheed Symposium on 
Magnetohydrodynamics; to be published. 

Bonp, J. 1957 Phys. Rev. 105, 1683. 

Keck, J. C., Kiver, B. & WENTINK, T. Jr. 1957 AICO Research Laboratory, Research 
Report no. 8. 


Lams, L. & Lin, S. C. 1957 F. Appl. Phys. 28, 754. 
PetscHEk, H. E. & Byron, S. 1957 Ann. of Phys. 1, 270. 


TURNER, E. B. 1956 Engineering Research Institute, University of Michigan, Rep. no. 
2189-2-T. 











195 


Slow viscous flow past a sphere in a cylindrical tube 


By HOWARD BRENNER and JOHN HAPPEL 


Department of Chemical Engineering, New York University 
(Received 29 November 1957) 


SUMMARY 


A theoretical treatment is presented for the slow flow of a 
viscous fluid through a cylindrical container within which a small 
spherical particle is confined. ‘The sphere is situated in an arbitrary 
position within the cylinder and moves at constant velocity parallel 
to the walls. Approximate expressions are derived which give 
the frictional drag, rotational couple, and permanent pressure drop 
caused by the presence of this obstacle in the original Poiseuillian 
field of flow. The primary parameters involved are the ratio of 
sphere to cylinder radius and fractional distance of the particle 
from the longitudinal axis of the cylinder. With appropriate 
modifications, the results are also applicable to a sphere settling 
ina quiescent fluid. ‘This yields the necessary boundary corrections 
to Stokes law arising in connection with devices such as the falling 
ball viscometer when the sphere is eccentrically located. 


1. INTRODUCTION 


The problem of viscous flow through a cylindrical duct containing 
particles of approximately spherical shape is of interest in connection with 
a variety of processes, such as fluidization, elutriation, and flow through 
fixed and moving beds of solids. A logical start towards elucidation of the 
hydrodynamic behaviour of these systems is undertaken here by considering 
the slow translation of a single spherical particle moving parallel to the 
longitudinal axis of an infinitely long circular cylinder through which a 
viscous fluid is flowing. The sphere may occupy any preassigned position. 
This is in contrast to the work of previous investigators who have concerned 
themselves exclusively with the axisymmetrical case in which the sphere is 
restricted to the cylinder axis. Haberman (1956, 1957) has recently reviewed 
the literature on this subject. 

In a companion article by the authors (Happel & Brenner 1957), the 
results of the present theoretical investigation have been employed to discuss 
the behaviour of dilute fluidized beds at low Reynolds numbers. 


2, NOMENCLATURE AND BOUNDARY CONDITIONS 
In the treatment to follow, it is necessary to resort to a variety of different 
coordinate systems. These are: Cartesian coordinates (x,y,z), spherical 
coordinates (r,4,¢), and cylindrical coordinates (p,¢,z), each having a 
N2 
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common origin at the sphere centre. It is also necessary to utilize Cartesian 
coordinates (.Y, Y, Z), and cylindrical coordinates (R, ®, Z), both originating 
ilong the cylinder axis and chosen so as to make s = Z. Relations between 
these various systems of coordinates are depicted in figure 1. 
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Figure 1. Relations between different coordinate systems emploved. 
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Figure 2. Definition sketch. 
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The sphere moves with an arbitrary constant velocity U relative to the 
cylinder wall in the direction of 3-positive, parallel to the cylinder axis. 
At the same time, the fluid flows in laminar flow with a mean velocity of } U4), 
in the same direction. ‘The sphere radius is a, the cylinder radius is Rp, 
and the centre of the sphere is situated at a distance 5 from the cylinder 
axis, as shown in figure 2. 

In terms of a coordinate system which moves with the sphere, the usual 
hypothesis of no relative motion at fluid-solid interfaces results in the 
following boundary conditions which define the fluid velocity field v: 


v=0 st 6r>@, (2.1) 

and 
v=-i,U at R= R,. (2.2) 
At large distances from the sphere, s = + &, the disturbance propagated 


by the sphere vanishes and the fluid velocity distribution becomes 
Poiseuillian. This gives rise to the additional boundary condition, 


v=1,{U,(1—R?/R?)—U} at z= + ow. (2.3) 
The equations of motion to be satisfied are 
V°v = (1/n)Vp, (2.4) 
together with the continuity equation for incompressible fluids, 
V.v=0. (2.5) 


Here yp is the viscosity and p the viscous pressure. Use of the linearized 

creeping motion equations restricts the validity of the final results to 

situations in which the relative particle Reynolds number, 
2a'U,(1—6?/R?)— U /n, 

is small, where 7 is the kinematic viscosity. 

As in a previous study (Happel & Byrne 1954), the above boundary- 
value problem is solved by the method of ‘reflections’. ‘Thus, the solution 
consists of the sum of a series of velocity fields, all of which satisfy (2.4) 
and (2.5), and each partially satisfying the boundary conditions as follows: 


v = i,{U,(1— R?/R2)— U}, (2.6) 
—y! at r= 4, n 

y) = . (2.7) 
G atse= to he r= co), 
-v) at R=R,, 

v2 = . (2.8) 
0 at Z= t. 

3 -y? atr=a 

Wr = 4 : (2.9) 

Q at 2 £ co (Ke: % = co}; 


etc., with as many fields taken as needed for an appropriate degree of 
approximation. ‘The field v satisfying the boundary conditions (2.1) to (2.3) 
is then obtained in the form 


V = VO+ yO 4 y+ yO84 (2.10) 
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and the corresponding pressure field is given by 


p= p+ pO + p® + pO +... (2.11) 
The success of this method of solution depends upon the linearity of the 
equations of motion. Its advantage resides in the fact that it is only necessary 
to consider boundary conditions associated with one surface at a time. 
If we let W and L, respectively, represent the frictional force and 
rotational moment experienced by the sphere, there are analogous relations 
of the form 


bo 


W = W® + W + W? + W® + ..., (2.12 a) 
and L = LO+L%+L0+L®+.... (2.12b) 


Furthermore, if AP, represents the additional pressure drop (above that 
due to the original Poiseuillian field, AP)) experienced by the fluid as a 
result of the presence of the sphere, then 

AP, = AP, +AP,+AP,+AP,+..., (2.13) 


where the pressure drop, AP,, associated with field 7 is simply the difference 


between the viscous pressure p” at Z= — © and Z= + ow. There is 
no ambiguity in this definition of pressure drop, since as Z+ + ~ the 


pressure becomes constant across the tube and, hence, independent of ® 
and R. It will prove convenient, therefore, in calculating pressure drop 
to evaluate the pressures p” at R= R,. Except for an arbitrary additive 
constant, which can be taken to be zero without loss of generality, the 
conditions of symmetry about the plane Z = 0 demand that the pressure 
be an odd function of Z. These observations lead to the formula 


AP; = —2 lim (p®)p =p,» (2.14) 
J+ 


where the subscript denotes evaluation of the function at the cylinder 
wall. 


3. "THE FIRST REFLECTED FIELD 
‘The unperturbed field, v, is defined in (2.6). It is easy to demonstrate 
that 
Ww) on 0, (3.1) 
and L® = 0. (3.2) 
The first reflected field v'”, defined in (2.7), has already been obtained by 
Simha (1936) in connection with suspension viscosity. Using Lamb’s 
(1932, p. 594) general solution in spherical coordinates, Simha’s result 
can be expressed in the form 
v) = V x (rx?,) + V(O®, + OF), + OY, + (1/2p)r?2p, — (1/30p.)r?2p™,) + 
+ (1p )e(p2 + $p",+ 4p), (3.3) 


and 


po = po +p, +p, oA) 


where r is the radius vector drawn from the sphere origin, and p_;,.4) 
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M_.,,.,) and y_¢,4,) are solid spherical harmonics of degree — (m+ 1) defined 
as follows: 
pw = Ar-*P, (cos @), p”,/u = BRyr- cos dP} (cos 8), 
p'?,|4 = CR2rP,(cos 4), ), = DR? rP, (cos 9), 3.5) 
Jed 
®”), = ER? r- cos dP} (cos 9), @®) = FRir—'P,(cos 6), | ( 
x2), = GRyr sin dP} (cos A), 
with the constants given by 
A = §a[U— U,{1 — (6?/R?)-3 ) 
| 
B= —'aU4(b/Ry)(a/Ro)’, C = —3aU,(a/Rp)', | 
D = }a(a/Ry)?[U— U,{1 — (6?/R2) — 3(a@?/R?)}], 


(a?/R?)}], 


E = — $aU,(b/Ry)(a/Ro)', F = —}aU,(a/R,)®, | 
G = aU,(b/R,)(a/R,)?. J 





The quantities P,, (cos) and P’’(cos @) are Legendre polynomials of order n, 
and associated Legendre polynomials (of the first kind) of order m and 
rank m, respectively. For reference, the following values appear in 
equation (3.5): 
P,(cos @) = cos@, P}(cos @) = sin @, 
P}(cos@) = 3sin@cos#, P3(cos@) = }(5 cos? @—3 cos @). | 


The drag and rotational moments experienced by the sphere can be 
obtained from the relations 


W® = —4rV(r*p"),), (3.8) 
and L® = —8rV(r?yx)). (3.9) 


These relations apply quite generally to Lamb’s solution. Noting from (3.5) 
that 3p"), = wArcos@ = wAz and introducing the value of A from (3.6) 
results in 


W® = —i, 6zpa[U — U,{1 — (62/R2) — 2(a?/R2)}]. (3.10) 


In a like manner, application of the relation rsindsin@ = x yields 
L® = —i, 82xpa?U,(b/Ry)(a/Ry)- (3.11) 


If figure 2 represents a meridian plane passing through both the sphere 
origin and cylinder axis, the tendency of the couple is to rotate the sphere 
in a clockwise direction. 

We shall limit ourselves in the subsequent development .» an 
approximate solution in which the ratio of sphere-to-cylinder radius a/R, 
is small. This covers most situations of practical interest. The terms 
in (3.5) are such that a final expression for the drag, correct to zeroth and 
first powers of a/R,, can ultimately be obtained by retaining only the p"), 
harmonic in the present solution. Thus, in place of the previous results, 
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we now resort to the approximate solution 


v) = (1/2n)V(r2p%,) + (1/2) rp, (3.12) 
and 
p> =p), (3.13) 
where 
pp. jp Hr-*cos@ H(z ry, (3.14) 
with the constant H defined by 
H = 3a(U— U,{1—(6?/R*)\]. (3.15) 
To the present degree of approximation, the drag is now given by 
W) = -1, 6mpa[U — U,{1— (67/R5)}], (3.16) 


while the rotational moment remains unaltered. 


4. ‘TRANSFORMATION TO CYLINDRICAL COORDINATES 

In order to compute the field v®, defined in (2.8), it is necessary to 

establish the form taken by v" at the cylinder wall, R= Ry. This is best 

done by transforming the latter field to cylindrical coordinates (R, ®, Z), 

originating along the cylinder axis. With the expression for Pp, u, given 
in (3.14), and the aid of the relation 


V(1/r) = -r/r’, (4.1) 
it is easy to show that an alternative form for the first reflected field is 
v) = LH[2(1/r)i. — V(z/r)]. (4.2) 
In terms of the desired coordinates we have 
i,= i, (4.3) 
and 
. © —% Be . 2 
VY =i, =p? - Rip * 1, =F (4.4) 


so that our objective may be attained by transforming the scalar functions 
1/r and z/r to these coordinates. 
Watson (1922) gives the relation 


1/r = (p? +27)—"2 = (2/77) | K,(Ap)cos Az da, (4.5) 


where K, is the modified Bessel function of the second kind of order zero. 
Although this transformation is only valid for p > 0, the restriction is 
inconsequential since it is sufficient for our immediate purposes to evaluate 
v only in the vicinity of the cylinder wall. Noting that 

p (R? + b?- 2hR cos D)!?, (4.6) 


we can avail ourselves of the further transformation, 


K,(Ap) = ¥ K;,(AR)I,(Ab) cos k®, (4.7) 


given by Watson. Here, J, and K,. are modified Bessel functions of the 
first and second kinds, respectively, of order k. The relation is valid only 
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for K>4b. For the reason outlined above, this restriction is of no 
consequence. ‘hese expressions combine to give 


Ibo 


I ‘nae 
-= > cosk® | K,(AR)I,(Ab)cos AZ di. (4.8) 
r 77 i “0 

‘lo obtain the second scalar transformation we observe, upon performing 
the indicated differentiations in the following equation, that the following 
identity is valid: 
sK,(AR)I,(Ab)cos Az = = [K,.(AR)I,(Ab)sin 1Z]—sin Mx [K,(AR)I,(Ab)]. 
0, Ci 


Thus, if (4.8) is multiplied by z, we obtain with the aid of the above 


~ p. 
~==5 SS cosk®[K,(AR)I,(Ab)sin AZ]! ~~ — 
r : 


YI bo 


J oO ra) ’ ; ; ; 

S cosk® | — [K,(AR)I,(Ab)]sinAZ da. 

—— ly OA 
h f ~ O 
The first term in brackets vanishes at the upper and lower limits of 
evaluation, A= (0 and x. Upon performing the indicated differentiation 
in the second expression, we arrive at the requisite transformation, 
~ ? of 
: =-- S cosk® | [RK/(AR)I)(Ab) +6K,(AR)L(Ab)]sin AZ da. (4.9) 
The differentiations denoted by primes are with respect to the entire 
argument. 

It is now a relatively simple matter to obtain the desired transformation, 


[vp ip(1/7) > cosk® | 4,(A)sin AZ da + 


-i,(1iz) ¥ sink® | 8,(A)sinAZ das 


+ i,(1/7) S cosk® | y,(A)cosAZ da, (4.10) 


where we have written 
a,(A) = H[AR, K,(AR,)L).(Ad) 4 
+ ADK (AR, )T,(Ab) + {R?/ (ARy) | K (AR, )L;,(A8)], (4.11) 
B,(A) = — Hk[(b/R,)K (AR, )L(Ab) + KAR )L;.(A5)], (4.12) 
y,(A) = HAR, K; (AR, )I;.(Ab) 4 
+ ADK (AR, Ti (Ab) + 2K) (AR) T;,(Ab)]. (4.13) 
These relations have been simplified to some extent by the use of Bessel’s 
modified equation, 
K/(AR) = — K,(AR)/(AR) + {1+ 4?/(Q?R?)} K, (AR). (4.14) 
An analogous expression for the initial pressure field at the cylinder wall 
can be obtained from (3.13), (3.14) and (4.8) with the assistance of the 


relation z/r3 = — e(1/r)/ées, yielding 


[p?], = (2uH/z) > cosk® | AK,(AR,)1,(Ab)sindZ da. (4.15) 
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5. THE SECOND REFLECTION, v‘~) 
A general solution of (2.4) and (2.5) in cylindrical coordinates, suitable 
for the field v®, is given by 


v= > | Vx G07) + HP +R V+, SF | (5.1) 
h L Ee 


is oR ‘ vA 
and 2 o2t]') 
a, es 5:2 


where Q,, ‘’,, and II, are arbitrary cylindrical harmonic functions of order k 
in (R,@, Z) coordinates. That this is, indeed, a simultaneous solution of 
the creeping motion and continuity equations is best verified by making 
the substitution Ro/OR = X(e/0X)+ Y(e/eY) in (5.1) and carrying out 
the operations indicated by (2.4) and (2.5) in Cartesian coordinates. 

These equations are perfectly general and would be applicable for the 
fields v, v®, etc. In the present instance, the harmonic functions are 
assumed to have the form 


Hi) = —(1/n)cosk® | An, (A)I,(AR)sin AZ dA, 
W2 = —(1/n)cosk® | A—Wp,(A)L,(AR)sin AZ dd, | (5.3) 
Q2) = —(1/m)sink® | A-tew,(A)I,(AR)sin AZ dA. 


U0 
The functions 7,(A), u%,(A) and w,(A) are to be determined from the 
boundary conditions expressed by (2.8). Upon substituting values from 
(5.3) into (5.1), we eventually find 


y2?) — —i,(1 77) > cos RD | | oe) 


A x 


kI,(AR) 


Co + ahy(A){(AR) + 





+a, (A NRIZOAR) sin a 


-ig(1/n) > sin kb a | 0) — m,(A)RI’ (AR) — 


fo NZ dd 


kI,(XR) 


4 





-i,(I/7) ¥ coskd | [w,(A)I,(AR) +7, (A)ARI.(AR) + 
+ m,(A)I,(AR)]cosAZ dd. (5.4) 


Evaluation of the above at R = Ry, and comparison to (4.10) by means of 


the boundary condition, 
2)]  — 1) = 
[v?], = —[v®],, (5.5) 
leads to three simultaneous equations involving the functions 7,(A), ;.(A) 
and w,(A) in terms of z,(A), 8,(A) and y,(A). Solving these equations and 
introducing values for the latter functions from (4.11) to (4.13) gives: 


7.(A) = HK,(ARp)1,(Ab){T,.(ARy)} 2 + E,(A), (5.6) 
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by(A) = — [1+ ARo (AR TARo)} 7 7A) + 
+ H{I,(ARy)}1[ARy K{(ARy)1 (AB) + 
+ ABK;.(ARg I; (Ab) + 2K;.(AR,)I;,(Ab)], (5.7) 
and 
w,(A) = 2R{AR, L(ARo)} tL (ARo)7(A) — HK; (AR )Z;.(Ab)], (5.8) 
where €;(A) in (5.6) is given by 
E.(A) = ALL, (Ab) T,(ARo ){L,.(ARo)} + — (6/ Ro )I.Ab)] x 
x [2k713(ARg){(ARo)*Z;.(ARo)} 1 + {ARy + (R?/(ARp)) 2 (ARy) — 
— 21;,(AR)L,(AR,) —ARo TAR»)? ]-1. (5.9) 
These equations have been simplified to some extent by means of the 
relations 
Ti (ARy) = —I,(ARo)/(ARo) + {1 + 2?/(A2R?)}T,.(AR»), 
and 
K,(ARo)I;.(ARo) — K;,(ARy)I,(ARo) a 1/(ARo). (5.10) 
The additional boundary condition, v®—>0 as Z+ + , is met because 
v +0 as Z> + o and v” and v® are related by (5.5). 

The pressure drop due to the first two reflected fields can now be 
computed. Substitution of the value of II, given in (5.3), into the 
expression for p”, given by (5.2), and evaluation at R= Rj, using the 
value of 7;.(A) from (5.6), yields 

16 mn <. on . or 
[P?]n= —-—— & cosk®| AK;,(ARo)I;,.(Ab)sin AZ dd— 
7T x /0 

aa 
= = cos k® | 


T } 


"NE, (A), (ARp)sin AZ dd. (5.11) 


( 
As evidenced by (4.15), the lead term in the above expression is simply 
—[p™],, so that 

[pP +p? Jp, = —(2u/7) SY cosk® | 2€E,(A),(ARy)sinAZ dv. (5.12) 
From (2.14), the pressure drop due to these first two reflections is, 
therefore, 


(AP, +AP,) = - a lim 3 cos k® [ none: 


0 Zen k Dn 0 A 





1A, (5.13) 
where we have written 
ny(A) = (ARy)PE (AT; (AR»)- (5.14) 
By Dirichlet’s theorem, 
— [ . (A) sinAZ A= 


/ 


nj (9+ ). (5.15) 





br! 


This limit, as A+0+, can easily be obtained by expanding the Bessel 
functions in series for small values of their arguments; whence, 
( —4H(1-—6?/R?) for k=0, 


0 for k # 0, 


m(0+) = 9,(0) = - (5.16) 
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from which, with the definition of H in (3.15), we finally obtain 
(AP, +AP,) = (12pa/R*)(1 — b?/ R2)[U,(1 — 6?/R*) — U). (5.17) 

A more physically meaningful representation 1s 
[((AP, + AP,)7R?2]}U, = [6zpa{ U,(1 —b?/R*) — U}|[U,(1 — 6?/R?)}. (5.18) 
For the small values of (a/R,) for which the present approximation is valid, 
the first term in brackets represents the additional force (above Poiseuille’s 
law) required to force fluid through the tube. ‘The second term 1s the 
mean velocity with which fluid traverses the cylinder. ‘Their product gives 
the additional energy dissipation incurred by the presence of the sphere in 
the original field of flow. On the right-hand side of this expression, the 
first term in brackets is the frictional force experienced by the sphere 
(equation (3.16)), whereas the second bracketed term corresponds to the 
local velocity of the unperturbed parabolic field in the vicinity of the sphere. 
It appears then that, for 2 sufficiently small sphere, the product of the drag 
and J/ocal velocity also gives the additional energy dissipation caused by the 
presence of an obstacle in the field of flow. 

Within the volume of space presently occupied by the particle, the 
field v'? has no singularities. ‘The sphere can therefore experience no 
frictional forces or couples in virtue of this field, and we have 


72) — () (5.19) 


’ 


and L® = 0. (5.20) 


6. "THE THIRD REFLECTION 

‘lhe results obtained thus far for drag and pressure drop may be regarded 
as the seroth approximation in powers of a/R,; that for rotational moment 
is correct to the first power. ‘To evaluate each of these quantities correctly 
to the next highest powers in a/R, it is necessary to consider the third and 
fourth reflections. Fortunately, exact solutions for v® and v™ are not 
required to arrive at exact values for these initial corrections. 

The frictional resistance and couple associated with v® can be calculated 
exactly by means of Faxen’s (1927) laws (see, also, Peres (1929)). In the 
present application, these laws take the form 

W® = 6zpa[v™], + 7a°[Vp Jo, (6.1) 
and L® = 47pa®[V xv], (6.2) 
where the subscript zero indicates that the function in brackets is to be 
evaluated at the sphere centre. ‘The relations are valid for arbitrary fields, 
v” and v, provided that each is in accord with the equations of slow 
motion and they are related by (2.9). ‘To preserve the consistency of the 
present results regarding powers of a/Ry, we omit the last term of (6.1). 

In terms of the (R, ®, Z) system of coordinates, the centre of the sphere 
is situated at the point (R = 6, ® = 0, Z=0). Thus, from (5.4) and (6.1) 


Ww i, 6u(a’R,) > | [ub (A)L,,.(Ab) 4 7, (A)AbT, (Ab) + 
| | a, (A)I,(Ab)]d(AR,). (6.3) 
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Employing the values of %,(A) and 7,(A), given by (5.6) and (5.7), and 
putting x = AR,, this becomes 





W® = —i, 6mpa[U — Uy(1 —02/R2)]f(b/Ry)(a/Ry), (6.4) 
where 
| ey. 3 & [fe oP @M@B) — arr ar]? 
} — b _ —— S § O7(% a BI iB Ba 
f(b Ry) rAY ) Qa — . > Erol I ,.(x) (28) 
eo) ‘a f 
. [ue vy] - 2Kj(%)Li(28) [1 ;.(%8) +%BI,(«B)] + dx, 
I ,.(%) . 


(6.5) 
in which 
Sts) = eee A (1 " wa )Fi(e) - 20 (%)Mi(%) _ UF | ' (6.6) 
a3 T;.(x) a2 o 
and for brevity we have put 8 = (6/R,). Equation (5.10) has been employed 
in simplifying the above. 
Using the method indicated in the Appendix, f(4/R,) has been evaluated 
as a power series in even powers of b/R, and the first two terms, obtained 





by numerical integration, are 
f(b/Ro) = 2:1047 — 0-6977(b/R,)? +.... (6.7) 
This development is valid for values of b/R,—0, near the cylinder axis. 
In the vicinity of the cylinder wall, where b/R,—>1, it is possible to make 
an exact calculation, thereby obtaining (see the Appendix) 
lim (1—6/R,)f(o/Ro) = 3. (6.8) 
b/ Rol 
This limit is approached asymptotically. A tentative plot of the function 
(1—5/R,)f(b/Ro) vs b/Ry, derived from these limiting expressions, is 
presented in figure 3. This should be of some assistance in the extrapolation 
of values of f(b/R,) beyond those for which (6.7) is strictly applicable. 
In view of the inadequacy of the formulae for intermediate values of the 
eccentricity b/R,, such a plot can only be regarded as preliminary. 
Equation (5.4) results in the expression 
[V xv], =i,(1/7) SY | [2Am,(A)i(Ab)—kb-1w,(A)T,(Ab)] dA, (6.9) 
k 0 (0) 
where we have noted that 


lig le = i,. (6.10) 
Substituting in (6.2) and employing both (5.6) and (5.8) eventually gives 
L® = —i, 87pa?[U — U,(1 —5?/R?)]2(b/R,)(a/Ry), (6.11) 


where 


g(0[Ro) = (8) = —(3/2n) z 1, pel aT (aB Mi (aB){Ta(2)) + 


+(e Kee) Serey IL Tite) 


“ Br, («8) |} dx. (6.12) 





The function 5,(x) is defined in (6.6). 
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Figure 3. Tentative plot of equation (6.5) for calculating drag and pressure drop. 
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Figure 4. Tentative plot of equation (6.12) for calculating rotational moment. 
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For small values of the eccentricity b/R,, numerical integration of the 
foregoing yields 
2(b/Ro) = 1:296(b/R,) +... , (6.13) 
and the expansion proceeds in odd powers of 6/Ry. In the limiting case, 
as the wall is approached, an exact calculation gives 


lim (1 —6/R,)?g(b/R,) = 2. (6.14) 
b/Ry>l 
Again, this limit is achieved asymptotically. A preliminary graph of 


(1 —6/R,)*g(b/R,) vs b/ Ry, prepared with the aid of these limiting forms, 
is shown in figure 4. 


7. ‘THE FOURTH REFLECTION 
In order to obtain the first correction to the zeroth approximation for 
pressure drop it is necessary to consider the field v. Here, again, an exact 
calculation to the first power in (a/R,)) can be made without a corresponding 
knowledge of v®. ‘This can be seen from the following argument. If we 
resort to Lamb’s general solution for the field v, reflected from the sphere, 
then by analogy to (3.8), we have 
4nV (rp?) = —W®, (7.1) 
This can be solved for the solid spherical harmonic function p™, by 
multiplying scalarly by r; whence, 


4cr eo (Fp) = - r. W®), 


or 
However, p'), is, among other things, a homogeneous polynomial of order 


—2. ‘Thus, by Euler’s theorem on homogeneous polynomials, 


This makes p>, = —(4nr*)*s. W. 
Using the expression for W® in (6.4), and noting that r.i. = z, yields the 
harmonic function, 
(1/n)p, = Jr-* cos 8, (7.2) 
where we have written 
J = 3a[U—U,(1 —6?/R?)]f (6/R,)(a/ Ro). (7.3) 
As in equations (3.12) and (3.13), that part of the entire velocity and pressure 
fields, v® and p®, associated with this particular harmonic function are 
v? = (1/(2n))V(r*p,) + (L/e rps, (7.4) 
and p = p™,. (7.5) 
Of all the terms in Lamb’s general solution, the only harmonic function 
which introduces the sphere radius a to the first power in the velocity and 
pressure fields is the p_, harmonic. ‘This is the dominant term as r—> «. 
Thus, as regards powers of a and, ultimately, in the next reflection, powers 
of (a/R,), this is the only term which need be retained to obtain results 
which are correct to the first power in (a/R)). To this degree of approxi- 
mation, (7.4) and (7.5) may be regarded as exact representations of the 
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velocity and pressure fields. Comparing these expressions to those given 
in (3.12) to (3.15), we find 


v3) = (J/H)v, p® = (J/H)p™. (7.6) 
lo the same degree of approximation, it is easy to see that 

v® = (J/H)v, p® = (J/H)p. (7.7) 
This immediately leads to the result 

(AP, +AP,) = (J/H)(AP, + AP,). (7.8) 


With the values for J and H, and the expression for the pressure drop due 

to the first two reflected fields, quoted in (5.17), this becomes 

(AP, +AP,) = (12a/R*)(1 — 6?/R*)[U,(1 — 6?/R?) — U) f(6/R,)(a/R,). (7.9) 
The fourth reflected field can have no singularities within the volume 

occupied by the sphere, so that 


WwW» =0, and L® = 0. (7.10) 


8. FINAL RESULTS 
Upon summing the individual results for drag, rotational moment, and 


pressure drop we find 


W i 6zpalU,(1 6?/R*)—U \{1 f(b/R,)(a/R,) + om (8.1) 
L = i, 8zpa*| U,(b/R,)(a/ Ry) - rl ‘(1 — b?/ R*) — U)g(b/R,)(a/Ro)? +... 
(8.2) 
and 
AP, = (12ya/R*)(1 — 6?/R*)(U,(1 — 6?/R=) — Ul[1+f(6/R,)(@a/Ro) +... ]. 
(8.3) 


These expressions are correct to the highest powers of (a/R,) quoted. 
The case of a sedimenting sphere in a quiescent fluid is obtained by putting 
U’, = 0 in the above. 

A form of the above equations suitable for examining situations in which 
the sphere is near the container walls can be obtained by expressing the 
previous results in terms of the ratio of sphere radius to minimum distance 
(of the sphere centre) from the wall, a/(R)—6). This results in 


W = i, 6zpa[U,(1 —52/R2) — uy +(1—b/R,)f( Ri) (3 “| +f 
(8.4) 








L = ~i, 8zpa? + U,(b/Ry)(1—5 Ro a a 
Ro- 





: ii. un i a ee | (9 5 
~[U,(1 —B2/R2) — U](1—8/R,)2e(b Ra) a5) Peeps (85) 


> = (12na/R2)(1 — b2/R2)[U,(1 —2/R2) — U] x 





: E +(1—5/R,)f(b RZ — ;) to | (8.6) 
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As previously observed, the quantities (1 — b/R,)f(b/ Ry) and (1 — b/R,)?g(b/ Ro) 
attain limiting values of } and 3, respectively, as b/R,->1. It should be 
noted that 
a 
Ro = 4, (8.7) 
where the equality sign applies when the sphere touches the wall of the 
cylinder. 


9. DISCUSSION OF RESULTS 

When the sphere is at the cylinder axis, b/R) = 0, equations (8.1) to (8.3) 
reduce to results previously given by Happel & Byrne (1954) and, 
independently, by Wakiya (1953) and Haberman (1957). 

It appears from the results of this investigation that both the drag and 
pressure drop are even functions of the eccentricity b/R). ‘This is to be 
expected, since symmetry requires that if the sphere is moved from its 
present location, R = 5, to the opposite side of the cylinder axis, R = —4, 
neither the direction nor the magnitude of the drag should be altered. 
Likewise, the pressure drop should be unaffected by this transition, in 
accord with present calculations. On the other hand, equation (8.2) shows 
that the rotational moment is an odd function of 6/R,. This conclusion, 
again, coincides with intuition, since the sphere will tend to rotate in a 
direction opposite to its original direction, without alteration in the 
magnitude of the rotational couple, when it is placed in mirror-image 
position on the opposite side of the cylinder axis. Of course, at the 
cylinder axis, we have 


L(0) = 0. (9.1) 


Further support for the validity of the present results is provided by 
observing that, in the absence of external forces acting on the sphere, the 
viscous pressure drop, due to the presence of a particle in the original field 
of flow, is an essentially positive quantity. ‘This stems from a one-to-one 
correspondence between pressure drop in rectilinear flow and energy 
dissipation. Thus, equation (8.3) shows that, for net flow in the 
+ z-direction, a positive pressure drop implies the inequality 


U,(1—b2/R2) > U. 


Inasmuch as U,(1—5?/R?) is the local fluid velocity, this shows that the 
sphere necessarily lags the fluid, an obviously correct inference. Moreover, 
equation (8.1) correctly indicates that, in these circumstances, the frictional 
force on the sphere is always in the direction of net flow, and the sphere 
is thus dragged along by the fluid. 

There are some remarkable conclusions to be drawn from the present 
calculations. The most interesting of these is, perhaps, the fact that the 
drag experienced by a small sphere sedimenting in a quiescent fluid does 
not increase monotonically as we proceed outward from the cylinder 
axis towards the wall. Rather, it attains a minimum value at some 


F.M. ° 
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intermediate point. ‘This can be seen quite clearly from (8.1) which, in 
present circumstances, has the form 

Ww i. 6mpaU[1 +f(b/R,)(a/Ry) +...]. (9.2) 
As we proceed away from the axis, the function, f(b/R,), and hence the drag, 
decreases initially as required by (6.7); however, as the wall is approached, 
equation (6.8) shows that f(//R,) ultimately increases like 1/(1—5/R,). 
Evidently, there is some intermediate value of eccentricity for which the 
drag is least. More extensive calculations are needed to locate this point 
with certainty. Experimental verification of the existence of a minimum 
may be a matter of considerable difficulty if it occurs near the axis, because 
of the relative magnitude of the ‘eccentricity’ term in (6.7) (for small 
values of (6/R,)) contrasted with the lead term. 

One might suspect that the results for the sedimentation of a small 
particle through a quiescent fluid, near the cylinder wall, would be exactly 
the same as for a sphere settling in the vicinity of a plane wall. Surprisingly, 
this is not the case. Lorentz’s (1907) results for the latter problem are 


: : 9 fa 
= —4.67 ; + —{—] +... 9.3 
Ww 1. 67pal I ia (5) } (9.3) 


while Faxen (1923) gives 


ae 2S. we. we 
L = —-1i,87pa?l E (;) _ sia (j) tn (9.4) 


where / is the perpendicular distance from the sphere centre to the wall. 
‘These should be compared with equations (8.4) and (8.5), respectively, 
which, when we set UU, = 0 and utilize (6.8) and (6.14), take the following 





form near the cylinder walls: 


W = -i. ompal| | + ia) ag | (9.5) 
ee a 0. 
and L i, 87a? Ll E (as) | (9.6) 


The distances from the wall, h and Ry—6, are comparable. Thus, in the 
case of frictional resistance, the difference between the two results is simply 
a numerical coefficient as regards the first power of the ratio of sphere 
radius to distance from wall. ‘The discrepancy in rotational moment is, 
however, of a more fundamental nature, since it involves differences in the 
exponents of the ratios, a/h and a/(Ry— 4). 

Experimental studies aimed at the verification of the foregoing equations 
are currently being conducted at New York University. 

The treatment of the present problem by means of the linearized 
Navier-Stokes equations fails to reveal the presence of sidewise forces 
tending to move the sphere towards the tube centre. In any real situation, 


‘Bernoulli’ forces tending to produce this result must exist. This lack 
of sidewise forces is a characteristic failing of the creeping motion equations, 
and stems from a neglect of fluid inertia in the original equations of motion. 
Since our present results show the force to be absent in the limiting case 
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as the Reynolds number tends to zero, it appears that, at low Reynolds 
numbers, this force should be proportional to the relative particle Reynolds 
number discussed in §2. Other considerations show this force to be an 
odd function of b/R,. 

It would be interesting to repeat the present study using the Oseen 
(1928) equations, which partially take account of fluid inertia. This, in 
conjunction with a corresponding experimental investigation of sidewise 
forces, would provide a stringent test of the Oseen equations. Here, the 
test would not be obscured by the fact that this approximate inertial effect 
is merely a second-order effect superimposed upon the primary viscous 
effect, as in Oseen’s (1928) correction of Stokes’s law. In this connection, 
it is of interest to note that the Oseen equations have already been employed 
by Faxen (1923) (see, also, Oseen (1928)), for the case of a sphere moving 
along the axis of a cylinder through a quiescent fluid. ‘This, of course, 
does not provide a test in the sense of the above. 
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APPENDIX 
Calculation of f(B) and g(8) as B->0 and B+>1 
The expression for {(8) given in (6.5) can be put in more tractable form. 
This is done by multiplying each side of the identity, 








Pin 9 
a [Ki(%)/1n(%)] = —1/[xlz(«)] (Al) 
by a/?(«8) and partially integrating the result, thereby obtaining 
” Ky(%) 7 | 2 eh 
I?(%B) + 2aBI,.(aB)I,(«8)] dx = | xlI?(«B) — + 
J, Fla | MOP AMP) ¢ aaa * 


x 
+ ["[ 4 n(%8) DY as 


The function in brackets vanishes at the upper and ues Sion and thus 


| 3° & (*(K,(a)I2(aB) 
i. 2 eee. 
MP=5 2 |. hee 








~ MFT (n)7 (xB) — BI la) Tp.q(aB) 2 bdo (A2) 
I? (x) J 


A series development of f($) can be carried out by expanding the Bessel 
functions of argument #8 in power se ries and evaluating the remaining 
integrals involving x by numerical techniques. An expression correct to 
powers in f* requires that we only retain the terms corresponding to 
k = 9,1 and —1 inthe infinitesummation. To this degree of approximation, 
we obtain 


F(B) = f(O) + $(B/7)(GQ1 — 202 + Qs — Qu). 
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The lead term, f(0), has already been evaluated by Happel & Byrne (1954) 
as well as by Haberman (1956), in connection with the motion of a sphere 
along the cylinder axis. ‘Their result (which has been checked) is 
3? [1—85(a)l?(«) o 
0) = —- ——————L—* | da = 2-1047. 
(9) | | P(x) % 7 


Numerical integration of the remaining integrals gives 


‘ iT) 


Qy= | PL Li(a)i* + Uta) ] da = 9-4488, 





Q=) = Ta) dx = —5-4094, 
x Pre 4 

on Bo (%)al (a) de mn os SOAS 

i—| a 2013, 


ej, Ae) 
‘These combine to give the expression for f(8) cited in (6.7). The modified 
Bessel functions of the second kind, which otherwise cause difficulty in 
numerical integration because of their behaviour at x = 0, have been 
eliminated by resorting to (A1) and variations thereof. 

To obtain the limiting form as B—1, we note from (A 2) that 
© Ky(a)IZGB) 
I;,(%) 

But, because of the monotonicity of the function J,(x) in conjunction with 
its essentially positive nature, one can demonstrate that 


° K,,(x)1?(xB) 


(1 —B)f(B) = 2 lim(1-B) > 


lin 
l k= oi) 


1 


K,{«a(1—B)]> > (A3) 
I(t FI k-—« I;,() 
by putting ® = 0 in (4.7) and (4.6). However, 

|. Ret Pp aa (A 4) 


from which the limiting form given in (6.8) follows without difficulty. 
A series expansion of g(f), defined in (6.12), correct to the first power 
of 6 gives 
: g(8) = — 4(B/7)(10, + 10s + 103), 
where 
ee ee 
~" Jy La) (1, («) + afe(«)] 
This results in (6.13). 
Finally, from (6.12), 





' - > iw 2? & (*% K,(a)l2z(aB) 
— 8)e(8) = ——lim(1-p??— 9S a s, 
ee eee ee 


Resorting to (A 3) and (A 4) it is now an easy matter to verify the correctness 
of (6.16). 
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The stability of a shear layer in an unbounded 
heterogeneous inviscid fluid 


By P. G. DRAZIN 
Trinty College, Cambridge 
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SUMMARY 
This paper considers the stabilizing effect of density 
stratification on the horizontal shear layer between two parallel 


streams of uniform velocities. A simple continuous velocity 
distribution, ( Vtanh(y/d), is used to represent the laminar 
shear layer. ‘The density of the fluid is assumed to vary as 


exp(— fv), v being the vertical coordinate, with a small total 
change in density across the shear layer. ‘The fluid is unbounded, 
and is assumed to be inviscid, incompressible and under the 
action of gravity. 

By the methods of hydrodynamic stability theory, it is shown 
that a disturbance of small amplitude and wave-number « is 
neutrally stable if the Richardson number, defined as J = g8d?/V?, 
has the value «?d?(1 — x?d?), and the form of the neutral disturbance 
is obtained. It follows that the critical Richardson number is }, 
so that the flow is stable if J > }. The relation between these 
results and Goldstein’s derivation of the same critical Richardson 
number for a flow with discontinuous velocity and density 
gradients is discussed. 


1. INTRODUCTION 

The natural occurrence of one parallel stream over another denser stream 
is widespread. Familiar examples are a river flowing into the sea and a 
warm layer of wind blowing over a cool one. ‘The knowledge that the 
tendency of gravity to keep the denser fluid below the lighter might strongly 
inhibit turbulence in such flows aroused interest in their stability many 
years ago. ‘This effect of gravity was first investigated by Stokes, who 
treated the stability of one fluid resting above another of greater density. 
Many authors extended Stokes’s work to various flows of several superposed 
horizontal layers of fluids with piece-wise constant density and horizontal 
velocity (such that the density and velocity are constant in each layer but 
may vary from layer to layer). ‘The stability of a heterogeneous fluid 
(i.e. a fluid with a continuous vertical variation of density) at rest also 
was treated by Lord Rayleigh and others. 

The general problem of stability of an inviscid fluid with both density 


and velocity continuously varying with height was approached in 1931 by 
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Taylor and Goldstein independently. ‘They each doubted the ability of 
a flow with several layers of homogeneous fluids to represent the phenomena 
occurring in the sea and atmosphere. ‘The mathematical question was 
whether the stability of a heterogeneous fluid could be approximated by 
use of alarge number of layers of homogeneous fluids. ‘To study this question 
they derived the equation of hydrodynamic stability of small amplitude 
waves in a parallel primary flow of an inviscid heterogeneous fluid under 
gravity. ‘They went on to consider the stability of some special flows with 
layers, in each of which either the shear (i.e. the velocity gradient) or the 
velocity is constant and the density either is constant or varies exponentially. 
They each concluded that a multi-layer system of homogeneous fluids 
could not be used to approximate the stability of a heterogeneous fluid. 

In § 2 of this paper the stability equation of ‘Taylor (1931) and Goldstein 
(1931) is derived in dimensionless form. In § 3 their results are summarized 
and discussed in the light of subsequent work. ‘The natural phenomenon 
of a shear layer in a heterogeneous fluid is more accurately described in § 4 
by use of continuously varying functions for velocity and density. This 
avoids the possible effect of discontinuities in either velocity or velocity 
gradient and either density or density gradient on the stability characteristics ; 
also the vorticity (i.e. the velocity gradient) 1s not assumed constant. 
It is then shown simply that a small-amplitude wave disturbance of 
(dimensionless) wave-number « is neutrally stable if the Richardson 
number J = «2(1—2?). It follows that the flow is stable for all wave 


disturbances if J > 4, the same critical Richardson number found by 
Goldstein (1931) for a shear layer with two discontinuities in the primary 


vorticity. 


2. ‘THE STABILITY EQUATION 

Following ‘Taylor (1931) and Goldstein (1931), we consider the inviscid 
stability equation for a parallel flow under gravity with a steady horizontal 
velocity U(y) varying with height y. ‘The fluid is heterogeneous with density 
p(y), although each material particle is incompressible. 

First consider the dimensionless parameters of this primary flow. 
Suppose that the velocity field has characteristic scales d of length and V of 
velocity ; and that the density distribution has length scale 1/8. ‘Then the 
flow may be typified by two dimensionless parameters. We define the 
Richardson number as 


J = gpd2/V?, (1) 


which measures the ratio of buoyancy forces to inertia forces. ‘The larger / 
the more stable the flow, because the steeper the density gradient the more 
energy of a disturbance is used to lift heavier fluid among lighter and push 
down lighter among heavier. We define the ratio of the length scale of 
velocity to that of density as 


L = Ba. (2) 
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It can be seen from the equation of motion that L represents the effect of 
variation of inertia arising from heterogeneity. Note that L =./F, where 
the Froude number F = V2/gd. 

The vector form of Euler’s equations of motion for a gravitating inviscid 


fluid is 





Du | 

— = —-~Vp—gVy. 3 

Di eo ee (3) 
The condition of incompressibility and the equation of continuity lead to 

Dp/ Dt = 0, (4) 
and V.u=0. (5) 
Now consider the velocity field 
u= U(y)+u’, v=v, w=), 

the pressure p= -g| pdyt+p’, (6) 
and the density p= p(y)+p’, 


representing the superposition of a primary flow and a small amplitude 
wave-disturbance (denoted by primes). 

The disturbance is assumed to be two-dimensional (so that w = 0) on 
the basis of the work of Yih (1955). Yih showed essentially that, for all 
parallel flows of a viscous heterogeneous incompressible fluid, the stability 
characteristics of a three-dimensional wave-disturbance are the same as 
those of a two-dimensional one at higher Richardson number J and lower 
Reynolds number. ‘Thus only two-dimensional disturbances need be 
considered in the search for a sufficient criterion of stability, with 
minimum J and maximum Reynolds number. 

Any sufficiently small disturbance may be resolved into Fourier 
components which are dynamically independent, and so primed quantities 
will be taken as a typical component, being the real part of a complex quantity 
proportional to exp{ix(x—ct)!, where « is a positive wave-number and 
¢ = ¢,+ic; is a complex velocity. ‘This wave has a phase velocity c, in the 
x-direction and logarithmic rate of increase xc; in amplitude; thus c; = 0 
for neutral stabilitv. ‘Taylor (1931) and Goldstein (1931) eliminated wu’, p’ 
and p’ from the equations of motion to get the equation for v’. Here 
equivalently we use the stream function of the disturbance 


os’ = o(y)exp{za(x—ct)}, (7) 
to integrate the equation of continuity with 
ue = ay’ oy, v= — oy’ lox = — tab’. (8) 


Then elimination of p’ and p’ from the equations of motion leads to the 


equation 


(U—e’ ~ 08) US — 2 Ces -V B= 0, 
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where primes now and henceforth denote differentiations with respect to y. 
‘This becomes Rayleigh’s equation if the fluid is homogeneous, i.e. if p’ = 0. 
When the appropriate dimensional scales are divided out of each quantity, 
this equation takes the dimensionless form 


(U—c)($" — 024) — Ud “pike $+L(p'/8p){(U—c)¢’ - U'd} = 0. (9) 


Taylor and Goldstein neglected the terms in L, which come from the 
effect of heterogeneity on the inertial terms in the equations of motion. 
This simplifies the mathematics, and L is in fact small in cases of practical 
interest where the density variation across the whole shear layer is small. 
When J 4 0, the type of the singularities of equation (9) is not altered by 
putting L = 0. 

In the cases examined by ‘Taylor and Goldstein the density is piece-wise 
or varies exponentially, i.e. p’/p = —f, where § is zero or constant in each 
layer. ‘[his restriction does not appear to alter any essential feature of the 
stability of the real flow in a shear layer. Equation (9) is thus transformed to 

(U —c)(p" — a6) — U"6+ J 4/(U—c) = 0. (10) 
3. DiIscUSSION OF THE WORK OF ‘TAYLOR AND GOLDSTEIN 

‘Taylor (1931) and Goldstein (1931) solved equation (10) for two cases. 

For constant U there are two exponential solutions 


b = exp[ + {a2 —J(U—c)-2}19]. 


For linear U the solutions were found in simple terms of Bessel functions 
of order }(1—4/)!? and complex argument. ‘l'aylor found 


d >(U—c)tttva 4J)} 


as U->c, if U is linear near the critical plane where U=c. ‘Therefore ¢ 
decreases to zero as U->c if 0< J <4; if J > }, then the singularity 
differs in that ¢ oscillates infinitely rapidly as well as decreases in amplitude. 
With either singularity the horizontal velocity of the disturbance 
(proportional to ¢’) is infinite. 

Taylor’s consideration of several problems of three and four superposed 
streams, each of homogeneous fluid, indicated that there might be stability 
for J > } in the limiting case of a continuous density distribution. This 
indication was supported by observation of the behaviour of streamlines 
depending on the singularity of ¢ near the critical plane, yet was apparently 
contradicted by his results for a continuous density distribution. In spite 
of the difficulty of handling the Bessel functions, he solved the problem 
with linear velocity and exponential density above a rigid horizontal plane. 
Only neutral waves can exist if J = g8/(velocity gradient)* > }, and no waves 
at all can exist if 0 <J < }. 

The condition J = } does not represent critical stability in the usual 
sense that all waves are stable for each / > } and that there is at least one 
unstable wave for each J <1. Rather, his solution implies that the 
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assumption of Fourier resolution ot an arbitrary disturbance into separate 
wave components is invalid for 0 < J < } in the inviscid flow considered. 
In fact, the possibility of a zero or negative critical Richardson number 
in ‘laylor’s semi-bounded flow can be shown by comparison with plane 
Couette flow between two rigid horizontal planes. Wasow (1953) proved 
that this flow of a homogeneous viscous fluid was stable to small disturbances 
at sufficiently large Reynolds numbers. ‘This means that the critical 
Richardson number for the inviscid stability of bounded plane Couette 
flow is non-positive, because density increasing with height is needed to 
create instability. ‘hus there may be some stable disturbances in 
addition to those waves found by ‘laylor tor J > U. 

Goldstein (1931) was principally concerned with a_ three-layer 
How having U =u, and p=p,exp(—2fh) for y2h; U =u,yjh and 
pP=p.expi|—P(y +h); tor h>y>-h; and U=—u, and p=p, tor 
-h>y. He found that disturbances can be neutrally stable only if 
J = gPh*/u> < }, and theretore the flow is stable or unstable according 
as J >} or J <j} respectively. 

Comparison of their results for piece-wise constant density with those 
for continuous density led ‘Taylor and Goldstein each to conclude that not 
all the essential features of a continuous flow could be approximated in 
the limit by a multi-layer system. ‘his conclusion was based on their 
particular results for three and four layers and on the doubt whether a 
multi-layer flow could approximate the singularity of the continuous flow 
in the critical plane where U = c. 

In 1951 Scorer derived equation (9) from the equations tor an n-layer 
system by letting m-» x. Also, Benton (1953) found the equation for the 
eigenvalues of the complex wave-velocity of an n-layer system and then 
letn—> «. For certain velocity and density distributions he found that the 
progressive wave-speeds were the same in the limit as those of the 
corresponding continuous flow. Recent authors have generally favoured 
the conclusion that the eigenvalues and eigenfunctions of a general 
continuous flow are identical with the limits of the corresponding values 
and functions in an n-layer flow, but no proof has been given yet. In any 
event the limit cannot be uniform in ” and y near the critical plane in view 
of ‘Taylor’s observation that an n-layer system cannot have the infinite 
horizontal velocity associated with the singular behaviour of ¢. 


4. ‘THE STABILITY OF THE FLOW WITH SMOOTHLY-VARYING VELOCITY 
AND DENSITY 
Discontinuities of velocity and density gradients or of uniform shear, 
as in the flows considered by ‘Taylor (1931) and Goldstein (1931), are 
discrepancies from the real flow which may influence the calculation of 
stability. It might be thought that any continuously curved velocity 
profile would be more difficult to handle than a profile consisting of straight 


line segments, but this is not always so. With a profile shaped like tanh y 
the stability equation (10) can be transformed to a second-order differential 
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equation with four regular singularities. ‘he antisymmetry of the profile 
then enables the stability criterion to be found for each wave-number 
without using any detailed properties of the general solution of the equation. 
Following Curle’s (1956) work on the corresponding flow of a homo- 
geneous fluid, let 
U = tanhy (— 0<y< o) (11) 
in dimensionless form. (‘This profile is shown in figure 1(a).) ‘Thus the 
velocity scale V is half the velocity difference across the shear layer. ‘The 
length scale d is V divided by the velocity gradient at y = 0. Comparative 
studies of stability suggest that the exact shape of the velocity profile is 
not of much significance in stability criteria. However, the profile (11) 
can be well matched to that of a free boundary layer between parallel streams 
as calculated by Lock (1951), who supposed the streams meet in a line 
perpendicular to their flow and used Blasius’s equation for the flow 
downstream of the line. 


q 


Se ee 











(a) (b) 
Figure 1. (a) The velocity profile U = tanh y. (6) A typical profile of the density 
p = exp(—Ly). 


Let p = exp(—Ly) (- o<y< ow) (12) 
in dimensionless form. The infinity of p at y= — « and the zero at 
y = + lead to a wave exponentially damped at y= + © where the 
primary flow is uniform, provided « > }Z (as can be shown by examination 
of equation (9) in asymptotic form). ‘There is no singular behaviour 
caused by the neglect of L, so equation (10) will be used here. A flow with 
negligible L is essentially equivalent to a real flow whose stability is 
determined by the wave-development in the interfacial shear layer, not 
in the surrounding uniform streams where all disturbances rapidly die away. 
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‘he inviscid boundary conditions that the vertical velocity vanishes at 
infinity become ap>0 as y>+t 0. (13) 
‘These conditions and equation (10) incidentally ensure that the horizontal 
velocity also vanishes at infinity. 

Before solving the eigenvalue problem, we note its time and space 
symmetry. ‘The equation and the boundary conditions are real because 
the inviscid flow is time-reversible. Therefore, if c is an eigenvalue 
corresponding to the eigenfunction ¢ for given values of « and //, 
c ¢,—te; 1s an eigenvalue corresponding to the function ¢* for the 
same ~ and.J. A nominal reverse of the y-direction changes the sign of g 
and f, but not of their product in.J. Also it changes the sign of U, an odd 
function of y, but leaves the boundary conditions unchanged. ‘There thus 
is no difference in the positive and negative v-directions in the specification 
of the problem. ‘Therefore, if ¢ is an eigenvalue corresponding to $(¥v) 
for given « and J, ~—c is an eigenvalue corresponding to ¢(—/) for the 
same « and J. ‘This gives, for each wave travelling in one direction, 
another similar wave in the opposite direction. On the grounds of 
uniqueness it is natural to suppose that, for a neutral disturbance with c; = 0 
at any rate, these waves coincide and the wave speed c, = 0. This certainly 
is true for a neutral disturbance of the flow with uniform density, i.e. with 
J =(, because Tollmien (1935) proved that c= U where U” = 0 for 
general monotonic LU’. The particular velocity profile U = tanhy is 
monotonic and is zero at the only point where U” = 0, namely at y = 0. 
Therefore, in seeking a neutral wave with c; = 0 for a criterion of stability 
of the flow of a heterogeneous fluid, it will be assumed that 

¢ = @, (14) 
This will be shown to lead to an eigensolution. ‘The uniqueness of this 
solution will be supported in the first paragraph of § 5. 
‘Thus equation (10) becomes 
U($" — a2) — U"b+J4/U = 0, (15) 
1.e.  +(2sech?y—«?+J coth? y)d = 0. (16) 
It remains to find eigenvalues J (x?) of equation(16) and boundary conditions 
(13). 

It is convenient to use U as the independent variable. It can then be 

shown that 
. ; | a? J 

it Mi , —1-@ O- ws? 
and af(U)=0 at U= +1, (18) 
where the subscript U denotes differentiation with respect to U. Equation 
(17) has four regular singularities, at —1, 0, 1 and «. In fact, d can be 
represented by the Riemann symbol 

oe () —1 

4(1+A) Wr - i, (19) 
jp Sy af 


” 


= @, (17) 
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where 
A= (1-4), p= +(a2—J)*2, (20) 
unless J = 0, in which case the symbol is 
( —1 1 20 ] 
P| la Ja 2 U} 
| —la -—da -1 | 


Consider the special case J =0 first. In this case a solution of 
Rayleigh’s equation for a homogeneous fluid is sought. Equation (17) 
becomes the associated Legendre equation of degree one and order «. 
These Legendre functions are 

¢ = (U—a){(1+U)/(1—U);'* or (U+a){(1— U)/(1+ U)}* 
in general. When « = 0 or 1, these two solutions coincide; the second 


solutions are then 





$= {Ulog >> -1 (x = 0), 

1+U U | 
a = jayne 1 ak q _ 
@ = (1-— U?) 2 oS a iy (a = 1) 


It can be seen that there are only two solutions satisfying the boundary 
: ying ; 
conditions (18), namely 
a= 0, ¢ = U = tanhy, ] 
tn 7 (21) 
and a= I, @ = (1— U?)'? = sechy. 

A natural way to solve the problem for small J would be to perturb 
these eigensolutions by use of a power series in J or some other convenient 
parameter. Before this can be done the singularities of 6 must be removed. 
They may be divided out by defining a new dependent variable 


x = (Ut 1) eU (1), (22) 
Then 
‘—-1 0 1 XL ] 
x=P| 0 0 0 f+e-1A U | 
satisfies 





4 ; 4 a a 
Xuv (4 : £ 1 > )xe + Gre — i - . , =0. (23) 
The boundary conditions are that x is regular at U = 0, +1. 

Note that the coefficient of the last term of this equation is of the form 
A(U+1)-"(U-—1)—.. The general form of this term for a function with 
four regular singularities at a, b, c, 2 (one exponent at each of the finite 
singularities being zero) is (A4U+B)(U—a)(U—6b) -(U—c)“, where 
A is the product of the exponents at infinity and B is some constant. Now 
the equation for 4, and therefore for x, is even in U. ‘This is why B = 0 
in equation (23) and the last term takes its simple form, 
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It now appears that it is unnecessary to expand x as a power series in J 
because an exact solution is given by 
x = constant, ($+u—4$A)(—$+p—-4A) =0. 
Now « > 0, and J is real. ‘Therefore the real part of A = (1—4/)"? is less 
than unity in magnitude. But the sign of » = +(x%?—./)!? has been chosen 
so that » has a non-negative real part. Therefore $+p— 3A 4 0 in any 
real flow. ‘Therefore 
—4+p-—]A=0. (24) 


On substituting for A and w from equations (20), it can be shown that this 
becomes 

J = «7(1—2?). (25) 
‘This equation is plotted in figure 2. The eigenfunction is found, by putting 
x = const. and J = x?(1—~2?) in equation (22), to be 


¢ = (sechy)**(tanh y)!-. (26) 
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Figure 2. The curve of neutral stability J = %?(1 —~”). 


This eigensolution, (25) and (26), joins up the only two solutions for 
J =. For each value of x, equation (25) gives a value of / for a neutral 
disturbance. ‘The two roots of « coincide if J = } and are complex if 
J >t. 

It may be noted that the horizontal velocity of the disturbance (i.e. 
d' exp(ixx)) is infinite in the critical plane where y = 0. ‘This would be 
modified by inclusion of viscous terms as is met in the asymptotic theory 
of the Orr-Sommerfeld equation for large Reynolds number. However, 
the velocity flux is finite because ¢’ = O(1/y) as y+ 0. Also, the vertical 
velocity of the disturbance, iafexp(iax), is finite everywhere if J is 


non-negative, 
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5. DIsCUSSION OF THE STABILITY CRITERION 

The relation (25) for neutral stability in the (/,«)-plane is shown 
in figure 2. ‘This is one part of the neutral curve, but other possible parts 
and the division of the plane into stable and unstable regions must be 
considered in order to see whether a wave is stable or unstable (i.e. whether 
c, < 0 or ¢; > 0) for any given values of J and «. The division of the plane 
will be made possible by the knowledge of the general stabilizing effect 
of the buoyancy forces and of the stability characteristics of a homogeneous 
fluid. In $4 it was shown that if J = 0 and c; = 0, then c, = 0, and that the 
only two wave-numbers for neutral stability are «=0,1. From the 
asymptotic viscous theory of the Orr-Sommerfeld equation it follows that 
the flow is stable for values (J, «) on the «x-axis if « > 1 and unstable if 
% <1. It may be concluded that the division of the plane into two simply- 
connected regions, the stable and the unstable, is as shown in figure 2. 
The possibility of a closed stable region to the left or unstable region to 
the right of the neutral curve must be admitted, but opposes the way in 
which gravity damps the energy of a disturbance and stabilizes a 
heterogeneous fluid. 

On these grounds the critical Richardson number, the maximum value 
of J for which there can be instability, is }. ‘This is the same as the critical 
value found by Goldstein (1931). ‘This striking agreement, despite the 
different density and velocity profiles used, naturally raises the suggestion 
that it is not accidental. It might be thought, for instance, on noticing 
that /’ = y near the critical plane y = 0 in both representations of the velocity 
profile in a shear layer, that the velocity near the critical plane effectively 
determines the stability of the flow. ‘The value = } obtained by Goldstein 
(1931) comes from the order of a Bessel equation, the order being principally 
determined by the behaviour of U near the critical plane. However, this 
behaviour does not appear to be so important for the smooth profile 
U =tanhy because the effect of the other singularities of the stream 
function of the disturbance can be seen in the analysis of $4. It has also 
been observed in $3 that plane Couette flow must have a non-positive 
critical Richardson number, although the velocity has the same behaviour, 


U'~¥, near the critical plane. ‘This shows that the value J = | is 
determined by the unbounded shear flow as a whole rather than by its 
critical plane alone. ‘Therefore, agreement, approximate or exact, is a 


justification of Goldstein’s use of his three-layer flow to determine the 
stability criterion. ‘The exactness of the agreement on a simple fraction } 
appears to be accidental, and is presumably a consequence of the simplicity 
of the two flows being compared. 

‘The condition / > } for inviscid stability is also sufficient for viscous 
stability, provided that viscosity is solely a stabilizing effect in this flow 
of an unbounded heterogeneous fluid. ‘The effects of density variation 
and viscosity are primarily separate, so it is possible to use flows of a 
homogeneous fluid as a guide. Heisenberg’s criterion is known not to be 
valid for unbounded flows, but there has been little detailed calculation 
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of their stability. Lessen (1949) calculated only part of the curve of neutral 
stability of a free boundary layer and found that a decrease of Reynolds 
number R caused an increase of stability in all his results. 

Lessen (1949) did not explicitly calculate the velocity profile of the 
boundary layer from Blasius’s equation, but Lock (1951) has done so. 
This permits a comparison of Lessen’s results with those of $4 after the 
dimensionless units have been related. ‘The inviscid asymptote corresponds 
very closely to x = 1, in agreement with the larger root of J = «?(1—«?) 
for J = (0, which shows that the stability characteristics of an unbounded 
shear layer depend only weakly on the exact shape of the velocity profile. 
This gives the anticipated position of the ends of the curve of neutral 
stability in the (x, R)-plane as./ increases from 0 to 4. The upper branch is 
surmized to drop from « = 1 to x = 1/¥2; a lower branch (not calculated 
by Lessen) similarly rises from «= 0 to x= 1/y2. When J =}, the 
branches coincide at x = 1/\2, and the flow is stable for all values of « 
and R. ‘This indicates that J > } is a sufficient condition of stability for 
a viscous shear layer between parallel streams of different densities. 


The author is very grateful to Dr G. K. Batchelor for much helpful 
advice and guidance during this work. In this time the author has been 
Twisden Student of Trinity College and a Research Student of the 
Department of Scientific and Industrial Research. 
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CORRIGENDUM 


‘* Heat transfer from surfaces of non-uniform temperature ”, by D. B. 
SPALDING (fF. Fluid Mech. 4, 1958, 22). 

Pages 29 and 30. In equations (13) and (14), and on figure 4, the 
coefficient 0-1 should be replaced by 0-2, 
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